A0-A188  934  A  DISCRETIZATION  OF  THE  INTEQIML  EQUATION  FON  THE  TINE 
DEPENDENT  LINEARIZ.  .  (U>  DAYTON  UN1V  OH  RESEARCH  INST 
K  0  OUDERLEY  AUO  87  UOR-TR-87-26  AFMAL-TR-87-38C1 
UNCLASSIFIED  F3M15-8E-C-228Q  FAS  1/1 


1*25 


AD-A188  b34 


AFWAL-TR-87-3061 


m  F'LE  COP'1 


A  DISCRETIZATION  OF  THE 
INTEGRAL  EQUATION  FOR  THE  TIME 
DEPENDENT  LINEARIZED  SUBSONIC 
POTENTIAL  FLOW  OVER  A  WING 


Karl  G.  Guderley 

University  of  Dayton 
Research  Institute 
300  College  Park 
Dayton,  Ohio  45469 


August  1987 


Interim  Report  for  the  Period  May  1986  to  June  1987 


Approved  for  public  release;  distribution  unlimited. 


FLIGHT  DYNAMICS  LABORATORY 
AIR  FORCE  WRIGHT  AERONAUTICAL  LABORATORIES 
AIR  FORCE  SYSTEMS  COMMAND 
WRIGHT-PATTERSON  AIR  FORCE  BASE,  OHIO  45433-6553 


DTIC 

^ELECTE 
^  \  QEC  0  8 1987 

O- 


87  11  30  102 


REPORT  DOCUMENTATION  PAGE 


la.  REPORT  SECURITY  CLASSIFICATION 

UNCLASSIFIED 


2».  security  classification  authority 


2b.  declassification/downgrading  schedule 


1b.  restrictive  markings 


3.  distribution/availability  of  report 

Approved  for  public  release; 
Distribution  is  unlimited 


4.  PERFORMING  ORGANIZATION  fjEPORT  NUM8£R(S) 

UDR-TR-87-26  1 

I 


6a.  NAME  OF  PERFORMING  ORGANIZATION  6b.  OFFICE  SYMBOL 

University  of  Dayton  at  applicable) 

Research  Institute 


6c.  AOORESS  (City.  State  and  Z IP  Code ) 

300  College  Park 
Dayton,  OH  45469 


5.  MONITORING  ORGANIZATION  REPORT  N  UMBER  (SI 

AFWAL-TR-87-3061 


7a.  NAME  OF  MONITORING  ORGANIZATION 

Flight  Dynamics  Laboratory  (AFWAL/FIBRC) 
Air  Force  Wright  Aeronautical  Laboratories 


7b.  ADDRESS  (City.  State  and  ZIP  Code) 

Wright-Patterson  Air  Force  Base 
Ohio  45433-6553 


Ba.  NAME  OF  FUNOING/SPONSOHING 

organization  A  ip  Force  .Wrights 
Aeronautical  Laboratories/ rlBR 
Air  Force  Systems  Command 


8c  ADDRESS  (City.  State  and  ZIP  Code) 

WPAFB  OH  45433-6553 


8b.  OFFICE  SYMBOL 
(If  applicable ) 

AFWAL/FIBRC 


9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 


F3361 5-86-C-3200 


10.  SOURCE  OF  FUNDING  NOS. 


PROGRAM 
E  LEMENT  NO. 


PROJECT 

NO. 


11  TITLE  (Include  Security  Clarification)  A  LM  S C Te 1 1  ZS 1 1  011  Of 

the  Integral  Eguation  for  the  Time  Dependent 


i2.  personal  authorisi  Linearize 

K.  G.  GUDERLEY 


13a.  TYPE  OF  REPORT  13b.  TIME  COVERED 

Interim  from  Nay  86  TOJune  87 


WORK  UNIT 
NO. 


611D2F 


ubsomc  Potential  Flow  Over  a  Wing 


14  DATE  OF  REPORT  tYr ,,  Mo  .  Day)  15  PAGE  COUNT 

August  19B7  138 


COSATI  COOES 


GROUP 


01 


04 


18.  SUBJECT  TERMS  ( Continue  on  reverse  if  necetsary  and  identify  by  block  number f 


Time  dependent  linearized  potential  flow,  integral 
equation;  fundamental  solutions  unsteady  aerodynamics. 


19.  ABSTRACT  (Continue  on  reverte  if  neceuary  and  identify  by  block  number t 


In  a  previous  report  the  integral  equation  for  linearized  time  dependent  subsonic  potential 
flow  over  a  wing  has  been  derived.  The  report  describes  how  this  formulation  can  be  dis¬ 
cretized.  The  wing  and  the  wake  are  divided  into  triangles.  The  flow  is  described  by  the 
time  dependent  potential  at  the  corners  of  the  triangles;  in  other  words  the  program  gen¬ 
erates  for  each  corner  a  table  of  the  potential  versus  time.  The  interpolations  necessary 
for  the  integrations  between  the  corners  points  and  also  with  respect  to  time,  are  describee 
in  detail.  Initially  the  upwash  is  assumed  to  be  known  as  a  function  of  time  over  the 
planform.  In  practice  one  will  express  the  upwash  by  a  rather  limited  number  of  standard 
functions  defined  over  the  planform  multiplied  by  a  superposition  of  step  functions  in  time 
Then  it  suffices  to  determine  for  each  of  the  spatial  standard  functions  the  response  to  a 
step  function  in  time.  In  this  form  the  method  can  be  used  in  aeroelastic  problems  where 
the  time  dependence  of  the  deformations  is  not  always  known  in  advance.  In  such  problems 
overall  forces  and  momenta  rather  than  the  detailed  pressure  distributions  are  needed. 


20  OISTRI8UTION/AVAILABILITY  OF  ABSTRACT 
UNCLASSIFIEO/UNLIMlTED  C33  SAME  AS  RPT  H  OTIC  USERS  □ 


22a.  NAME  OF  RESPONSIBLE  INDIVIDUAL 

Charles  L .  Kel ler 


21  ABSTRACT  SECURITY  CLASSIFICATION 


UNCLASSIFIED 


22l>  TELEPHONF  NUMBER 
t  Include  \  rro  (  'ode  > 

(SI  3)  255-7384 


22c.  OFF  ICE  SYMBOL 


nrwAi  n  ibrc 


.  +tm-  lk*m  A' A  «  1 t  ***  j'J|  *'«  <1 


iLiL  u, •  Ala* aL< iL'it/ii .,ik*it.*«Wt«L'iLkiLVilftL%Llii%LSIk^LHLl^i;1 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE 


19.  Abstract  (cont) 


'This  is  discussed  in  some  detail.  After  the  initial  response  to  a  step  function  the  flow- 
field  is  rather  smooth  (although  still  time  dependent).  Therefore  it  can  be  described  by  a 
reduced  number  of  parameters  and  one  can  proceed  in  greater  time  steps.  The  procedure  by 
which  this  can  be  done  is  also  described. 


UNCLASSIFIED 


Sf  ‘.URIIV  n  ASSir  ICATION  of  tm»s  page 


y; 


PREFACE 


The  work  was  performed  during  the  period  May  1986  through  January  1987  under 
contract  F33615-86-3200  entitled  "A  Study  of  Linearized  Integral  Equation  for  Steady 
and  Oscillatory  Supersonic  Flow”  to  the  University  of  Dayton  for  the  Aeroelastic  Group, 
Analysis  and  Optimization  Branch,  Structures  Division,  Air  Force  Wright  Aeronautical 
Laboratories,  Air  Force  Systems  Command  under  Program  Element  61102F,  Project  No. 
2304,  Task  Nl,  Work  Unit  22.  Dr  Karl  K.  Guderley  of  the  University  of  Dayton  Re¬ 
search  Institute  was  the  Principal  Investigator.  Dr.  Charles  Keller,  AFWAL/FIBRC  was 
Program  Manager. 


The  author  would  like  to  express  his  appreciation  for  the  excellent  typing  done  by  the 
staff  of  the  University  of  Dayton  Research  Institute,  in  particular,  Mrs.  Louise  K.  Farren 
and  Ms.  Roxanne  Sant. 


b 

b 

b 


TABLE  OF  CONTENTS 


SECTION  PAGE 

I  INTRODUCTION  1 

II  COMPILATION  OF  THE  RESULTS  OF  REFERENCE  1  2 

III  FORMULATION  IN  AN  UNDISTORTED  COORDINATE  7 

SYSTEM 

IV  SOME  GENERAL  REMARKS  11 

V  THE  SOLUTION  PROCESS  DESCRIBED  IN  14 

GENERAL  TERMS 

VI  INDEXING  OF  CORNERS  AND  OR  AREAS  18 

VII  THE  INTEGRATION  PROCEDURE  22 

VIII  REPRESENTATION  OF  h  FOR  DISTANT  26 

TRIANGLES  AND  TRANSFORMATIONS  REQUIRED 

FOR  THE  INTEGRATION  PROCESS 

IX  APPROXIMATION  FOR  h  IN  THE  AREA  46 

SURROUNDING  A  CONTROL  POINT  AND 
TRANSFORMATIONS  REQUIRED  FOR  THE 
INTEGRATION 

X  THE  NUMBER  OF  PIVOTAL  POINTS  TO  BE  89 

USED  FOR  THE  INTEGRATION  WITHIN  ONE 

TRIANGLE 

XI  SPECIAL  CHOICE  OF  THE  TIME  DEPENDENCE  92 

OF  THE  UPWASH 

XII  TRANSITION  TO  SMALLER  VECTOR  h  100 

XIII  EXAMPLES  FOR  THE  TRANSITION  TO  A  SMALLER  108 
VECTOR  h 

XIV  THE  METHOD  IN  THE  FRAMEWORK  OF  AN  116 

AEROELASTIC  PROBLEM 

XV  EFFECT  OF  A  TRUNCATION  OF  THE  WAKE  121 

REFERENCES  123 


LIST  OF  ILLUSTRATIONS 

PAGE 


Planform  124 

Modified  Planform,  Subdivision  into  Triangles  124 

Numbering  of  Triangles  125 

Numbering  of  the  Vertices  at  and  Adjacent  to  125 

the  Control  Point  ij 

Numbering  of  the  Corners  at  and  Adjacent  to  126 


a  Control  Point  on  the  Axis  of  Symmetry 
(j  -  0) 


Triangle  m  «  0,  n,  1  (at  the  Leading  Edge)  126 

Sweep  Angle,  l  n  and  uv  Systems 

Same  Triangle  as  Figure  6,  uv  and  uv  Systems  127 

Triangle  m  -  0,  n,  2  at  the  Leading  Edge,  127 

£  n,  uv  and  uv  Systems 

Triangles  m  -  i t - 1  ,  n,  1  and  m  =  i t - 1 ,  n,  2,  128 

Triangle  in  the  Wake  and  uv  Systems 

Triangle  n  =  i —  1  ,  n,  1,  uv  and  uv  Systems  128 

Triangle  m  -  i  -1 ,  n,  2,  uv  and  uv  Systems  129 

Vicinity  of  a  Control  Point  at  the  Trailing  129 

Edge  (i  -  i t )  i  «  it,  and  uv  Systems 

Vicinity  of  a  Control  Point  in  the  Interior  130 

at  the  Axis  of  Symmetry 

Control  Point  at  the  Axis  of  Symmetry  next  130 

to  the  Leading  Edge 


SECTION  I 


INTRODUCTION 

The  equation  for  three  dimensional  subsonic  unsteady 
potential  flow  possesses  fundamental  solutions  in  a  closed  form. 
Therefore,  one  can  readily  bring  the  problem  of  the  flow  over  a 
wing  into  the  form  of  an  integral  equation  with  two  independent 
variables.  This  formulation  is,  however,  not  directly  suitable 
for  numerical  evaluation,  for  it  requires  that  one  carry  out  a 
limiting  process  in  which  one  approachs  the  planform  from  above 
or  below.  In  Reference  1  transformations  have  been  carried  out 
which  generate  a  form  in  which  all  expressions  have  a  direct 
meaning  within  the  planform  of  the  wing. 

In  these  derivations  no  simplifications  (except  for  the 
initial  assumption  of  a  linearized  flow)  have  been  made.  For 
numerical  purposes  a  discretization  must  be  carried  out.  The 
present  report  describes  a  possible  procedure  of  this  kind.  Most 
of  the  report  deals  with  the  discretization  process.  A 
preliminary  code  for  the  two-dimensional  case  is  found  in 
Reference  2. 

The  method  proposed  here  is  rather  versatile  but  time 
consuming.  In  order  to  economize  one  must  know  which  information 
is  required  in  the  setting  in  which  the  results  will  be  used. 

The  author  has,  therefore,  added  some  sections  which  deal  in 
general  terms  with  the  form  in  which  these  results  will  be  used 
in  aeroelas t i ci ty .  This  includes  an  attempt  to  characterize  the 
necessary  time  resolution,  and  a  method  of  reducing  the  size  of 
the  problem  after  initial  "irregularities"  have  died  out. 


1 
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SECTION  II 


COMPILATION  OF  THE  RESULTS  OF  REFERENCE  1 

The  results  obtained  in  Reference  1  are  repeated  here, 
arranged  in  a  form  better  suited  for  the  numerical  work.  In 
Reference  1  the  Prandtl-Glauert  coordinate  distortion  has  been 
introduced.  For  unsteady  flows  it  is  less  useful  than  for  steady 
flows  because  it  obscures  to  some  extent  the  manner  in  which 
perturbations  propagate;  moreover  the  effect  of  sweep-back  is 
more  clearly  expressed  without  this  distortion.  We  shall  go  to 
the  undistorted  (but  dimensionless)  variables  in  Section  III.  To 
prepare  for  it  we  characterize  variables  in  the  distorted  system 
(which  in  Reference  1  appear  without  tilde)  by  a  tilde.  This 
will  allow  us  to  use  variables  without  tilde  through  the 
undistorted  system. 

In  this  section  we  compile  the  results  in  the  distorted 
system  of  coordinates.  In  Reference  1  we  denoted  by 

U  =  free-stream  velocity 

a  =  free-stream  velocity  of  sound 

M  *  U/a  =  free-stream  Mach  number 

p  =  free-stream  density 

L  =  a  characteristic  length 

The  coordinates  in  physical  space  and  the  time  are  denoted 
by  x»  y .  z»  and  t.  The  x  direction  is  the  free-stream  direction, 
the  planform  of  the  wing  lies  in  the  x,y-plane.  Let  <f(x,y,z,t) 
be  the  potential  and  p  the  perturbation  pressure.  Dimensionless 
quantities  and  the  Prandtl-Glauert  distortion  are  introduced  by 

x  =  x/L 

;■>  i  /•>- 

y  =  (1  •  M  )  "  y /L 

?  i  /  ’ 

z  =  (1  -  M  )  '  z/L 


V -.• »/ 


t  =  (at)/L 


--  P  -  1  /?  -  9  _  1  /  ?  _ 

♦  (x,y,z,t)  -  UL(j>(xL,(1  -  M  )  yL,(1  -  M  )  zL.tL/a) 

p  -  p/pu2 

Let  x  =  x  (y)  be  the  equation  of  the  trailing  edge.  Then 

x  =  xtr(y) 


with 


xtp(y)  -  L-1xtr((1  -  M^)"1/2yL) 

Let  z  =  g(x,y,t)  give  the  time  history  of  the  displacement  or 
deformations  of  the  wing  surface,  and 

-  ~  ~  ~  _i_  -  p  _i/p_  - 

g(x,y,t)  =  L  g(xL,(1  -  M  )  yL.tLa  ) 

To  express  the  boundary  conditions  at  the  wing,  Reference  1 
introduces 


w0(x,y,t)  »  (1  -  M2)"1/2(g_  +  (1) 

x  t 

The  dependent  variable  in  the  problem  is  h(x,y,t).  One  has 

$ ( x ,y ,z=±0, t )  =  ±2nh ( t , x  , y  ) )  (2) 

The  notation  z  -  ±0  means  that  one  approaches  z  =  0  from  positive 
or  negative  values.  The  values  of  h  in  the  wake  are  expressed  by 
its  values  at  the  trailing  edge.  To  emphasize  this  we  introduce 


Then 


h  ( t ,  x  ,  y )  =  htr(t  -  M  \x  -  xtr),y))  (3) 

In  Reference  1  the  notation  h^^,  h^2^  ,  and  h^)  has  been 
introduced  for  the  derivative  of  h  with  respect  to  its  first, 
second,  and  third  argument. 

In  the  integro-dif f erential  equation  for  h,  x  and  y  are 
points  of  the  wing  areas,  the  umbral  variables  \  and  n  range  over 
the  wing  and  the  wake.  For  each  choice  of  the  point  x,y  the 
integration  region  is  divided  into  an  area  surrounding  the  point 
x,y  denoted  by  A^,  and  the  remaining  area  (wing  area  +  wake  -  A.) 
denoted  by  A2«  The  contour  of  A^  is  denoted  by  .  One  has  the 
following  auxiliary  expression 

p  =  [(;  -  x)2  ♦  (n  -  y)2]1/2 

ret  =  (1  -  M<d~1[MU  -  x)  +  p] 

t  =  t  -  ret 

j(2)  _  _ 

p[M  U  -  x)  +  p]2 

j(3)  =  (1  -  M2 ) ( n  -  y) 

p[M  U  -  x)  +  p]2 

Then  one  finds  in  Reference  1  the  following  i ntegro- d i f f eren t i al 
equation  for  h. 

wQ(x,y,t)  =  -2n(1  -  M2)"1/2[h(1)(t,x,y)  ♦  Mh ( 2 } ( t , x , y ) ] 

+  B?  +  Bg  +  §9  (4) 


where 


B.j  and  By  are  obviously  closely  related  to  each  other,  but 
does  not  converge  in  a  region  containing  the  point  x  .  y  .  For 
details  see  Reference  1. 

If  h(t,x,y)  has  been  determined  up  to  a  time  t  for  all 
points  within  the  planfortn  (and  because  of  Eq .  (3)  also  within 
the  wake),  then  Eq .  (4)  allows  one  to  evaluate  9h/9t  = 
h^1  (t,x,y)  and  then  proceed  by  another  time  step.  After  h  has 

been  found,  one  evaluates  4>(x,y,z  =+0,t),  Eq .  (2),  and  the 
dimensionless  deviation  from  the  free-stream  pressure 

p  =  -(M*1*.  +  1)  _ ) 
t  x 

and  the  actual  deviation  from  the  free-stream  pressure 

P  =  pU"p 

The  factor  M  ^  in  the  formula  for  p  occurs  because  the  time  t  and 
4>  have  been  made  dimensionless  with  the  velocity  of  sound  a  and 
the  free-stream  velocity  U,  respectively.  (If  one  makes  t 
dimensionless  with  a,  then  a  perturbs  Mon  travels  by  the 
characteristic  length  L  during  a  time  interval  1  measured  in  t. 


If  one  makes  t  dimensionless  with  U,  then  a  particle  moves  by  the 
characteristic  length  L  during  the  interval  1  measured  in  t. 

Both  definitions  have  a  sound  physical  meaning.) 


5v»Vr 


SECTION  III 


FORMULATION  IN  AN  UNDISTORTED  COORDINATE  SYSTEM 


One  returns  to  an  undistorted  coordinate  system  (x,y,z,t) 
by  the  transformation 


x  =  x  /  L  =  x 


,  .  ..2,1/2-,.  ,  .  ..2,1/2 

y  =  ( 1  -  M  )  y/L  =  ( 1  -  M  )  y 


*.2,1/2-,.  ..2,1/2 

Z  =  ( 1  -  M  )  z/L  =  ( 1  -  M  )  z 


t  =  at/L  =  t 


5  =  C 


,  ,  2,1/2 
n  =  ( 1  -  M  )  n 


2  1/?  pi/"? 

<i>( x ,y , z , t )  =  0(x,(1  -  M  )  y ,  ( 1  -  M  )  "z,t) 


The  basic  data  in  physical  space  are,  of  course,  the  same. 
The  time  history  of  the  displacements  or  deformation  is  given  by 


z  =  g(x,y,t) 


Hence,  in  undistorted  coordinates 


z  =  g ( x , y , t )  =  L  g( Lx , Ly , Lt/a ) 


We  introduce 


htr ( t ’ y  )  =  h(t  ,xtp (y  )  ,y  ) 

Then  one  has  within  the  wake 

h ( t ,  x  ,  y  )  =  htr(t  -  M  ( x  -  xtr(y),y)) 
One  needs  the  following  auxiliary  functions 

r  ,  r  ,2  ,  ,  2 «  ,  >2 ,1 /2 

p  =  L  (  K  ~  x )  +(i-M)(n-y)] 


3 


rv*  wvawv  jr*.  vt”,y^\t"  t r»  ir»  ir 


,2.-1 


ret  =  (1  -  M  )  [M(C  -  x)  +  p] 


t  =  t  -  ret 


(2)  £  -  x  «■  Mp 


p[M(  £  -  x )  +  p] 


2 


.(3)  ( 1  -  m  )  ( n  -  y) 


p  [  M  (  £  -  x)  +  p  ]  ^ 


,(  3) 


The  expression  f  differs  from  the  expression  resulting  from 


substitution  into  f^3^  by  a  factor  (1  -  M2)^2.  After  multi- 


?  1  /2 

plication  by  ( 1  -  M  )  ,  one  then  obtains 


wo(x,y,t)  =  -2  tt  [  h  ^  1  ^  (t  ,x  ,y  )  +  Mh^(t.x.y)] 


*  B7  *  B8  +  B9 


(5) 


where 


B?  =  (1  -  M  ) 


(  2 ) .  (  2 ) 


( 3) u ( 3) 


f  h  (t , £ , n)  -  r J/hv (t ,£,n)d£dn 


(6) 


(1  -  «3)h<2,(t,x.y)  vk  -  <n 

J  o[M(F,  -x)+p] 

1  1 


h(3)(t 


,  x  ,  y  )  vj*  (  £  -  x  ♦  Mp  ) 

r. 


(_£  -  x  »  Mp)d£  +  (1  -  M‘~)(n  -  y)dr 
p[M  (  £  -  x)  +  p]2 


B_  -  ( 1  -  M  ) 


{f(2) (h(2) ( t , £  ,  n)  -  h(2) (t.x.y)  ) 


♦  f (  3)  (h  (  3)  ( t  ,£ ,  n)  -  h(3) (t.x.y)  )  )d£dn 


After  h(t,x,y)  has  been  found,  one  determines 


<f>(t,x,y,z  =  +0  ,  t )  =  2irh(t,x,y) 


p  =  -(M  <))t  +  <J>X) 


,.2  ,  x  y  ta. 
p  ( x  ,  y  ,  t )  =  pU  p(-j-,  y,  jj-) 


The  problem  will  be  discussed  in  the  form  shown  here.  In  the 
computation  we  evaluate  for  each  time  step  h^^  =  The  equation 

is  therefore  used  in  the  form 


3 1 ( t >x-y)  = 


-(2  it)  wq  (  x  ,  y  ,  t )  +  Mhx(t,x,y) 


(  2tt  ) _  [  B?  +  Bg  +  Bg] 


»rwewvwv%nm  ■  *t  i"* '  vmrv'  vrur  nr  v.  v  nr  ^  n.-  ^  v  ^  V>  W^A-  wuim  m '.  wjfuwv  m'jm 


SECTION  IV 

SOME  GENERAL  REMARKS 


It  is  desirable  to  hold  the  number  of  parameters  by  which 
the  wing  potential  is  described  small.  In  the  present  approach, 
this  means  that  the  number  of  points  within  the  planform  for 
which  the  potential  is  determined  should  not  be  excessive.  A 
lower  limit  for  the  number  of  points  is,  however,  unavoidable. 

By  discretizing  in  space  one  also  restricts  the  resolution  in 
time,  because  the  method  is  unable  to  follow  the  propagation  of 
pressure  waves  through  the  individual  surface  elements.  As  a 
consequence,  the  high  frequency  components  of  the  motion  are 
lost,  more  so  in  a  coarser  than  in  a  finer  grid.  A  minimum 
number  of  subdivisions  of  the  surface  is,  therefore,  necessary. 

At  the  leading  edge,  the  potential  shows  a  singular 

1  /2 

behavior.  In  a  steady  flow  it  behaves  as  u  multiplied  by  a 

constant  where  u  is  the  distance  from  the  leading  edge.  The 

-1/2 

pressures  then  behave  as  const  times  u  .  If  one  approximates 
the  potential  at  the  leading  edge  in  a  simple  manner,  namely  by  a 
linear  function,  then  the  pressure  in  this  region  will  be 
constant.  This  by  itself  is  not  objectionable;  in  the  interior 
of  the  wing  such  an  approximation  is  sufficient.  But  experience 
has  shown  that  the  integrated  pressures  (the  contributions  to  the 
total  lift  and  moment)  for  the  mesh  elements  at  the  leading  edge 
are  very  inaccurate.  One  might  remedy  this  by  a  correction 
factor  (universal  for  all  elements  at  the  loading  edge  and 
restricted  to  this  area).  This  possibility  might  be  worth 
exploring,  it  might  lead  to  a  simplification  of  the  programming 
work.  In  the  present  report  the  author  has  chosen  to  approximate 
the  potential  in  the  vicinity  of  the  leading  edge  by  a  function 
which  behaves  as  u1^  and  thus  to  anticipate  the  behavior  for  the 
steady  case.  Details  will  appear  later. 

In  the  present  report  the  analysis  is  restricted  to  the 
wing.  This  should  be  taken  as  the  first  step  towards  an  approach 
in  which  wing  and  fuselage  are  taken  into  account  simultaneously. 


For  a  slender  fuselage  a  linearized  approach  is  again  applicable. 
The  treatment  of  the  fuselage  is  simpler  than  the  treatment  of 
the  wing,  at  least  in  principle,  for  the  potential  can  be 
expressed  by  a  distribution  of  sources  over  its  surface  rather 
than  by  doublets,  and  the  singularities  in  the  kernel  of  the 
integral  equation  are  less  severe.  Tiieor  et  i  cal  1  y  ,  one  could 
apply  such  a  treatment  also  to  the  wing.  Then  one  would  not 
reduce  it  to  its  planform  but  consider  it  as  three-dimensional 
and  satisfy  the  boundary  conditions  separately  at  the  upper  and 
lower  surfaces.  The  treatment  in  terms  of  doublets  is,  however, 
preferable;  the  upper  and  lower  wing  surfaces  are  close  to  each 
other.  As  a  consequence  of  the  strong  interference  between 
sources  at  the  two  surfaces  the  result  will  appear  as  the 
difference  between  large  numbers.  This  difficulty  vanishes  if 
one  uses  doublets  rather  than  single  poles  at  the  wing  surface 
with  opposite  signs.  A  method  which  uses  single  poles  at  the 
fuselage  and  at  the  faired-in  transition  from  the  fuselage  to  the 
wing  and  doublets  at  the  planform  can  probably  be  developed 
without  difficulties.  In  each  representation  one  has  to  express 
the  velocity  component  normal  to  the  planform  and  to  the  fuselage 
in  terms  of  the  unknown  parameters  describing  the  doublet  and 
source  strengths.  From  these  data  the  potential  can  be 
computed . 


If  one  disregards  the  fuselage  one  is  inclined  to  take  the 
symmetry  of  the  problem  into  account  by  considering  a  planform  as 
shown  in  Figure  1.  This,  however,  generates  a  difficulty  at 


point  A  if  the  wing  has  sweep.  TSn-re  the  flow  pattern  becomes 


very  complicated.  The  u'  '  singularity  of  the  potential  at  the 
leading  edge  still  persists,  but  nu  analytical  solution  for  tho 
interior  is  available  which  could  be  used  in  the  computation  to 
anticipate  the  local  character  of  the  flow  field.  To  avoid  this 
difficulty,  we  shall  cons i dor  i  planform  is  shown  in  Figure  2. 


In  the  vicinity  of  the  wing  root  an  ipproaeh  which  replaces 
a  wing  by  its  planform  is  uns  it.  i  <>f  aelory  in  any  u-e, 
fairing  in  between  wing  a n d  fur-«i  i,'"  not  t  ikon  into 


The 

iccount  it 


all  and  the  effect  of  the  fuselage  on  the  wing  is  modeled  only 
approximately.  Therefore,  it  does  not  really  matter,  whether  one 
treats  the  problem  of  Figure  1  or  Figure  2.  A  satisfactory 
treatment  will  always  require  a  theory  in  which  the  fuselage  is 
included . 
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SECTION  V 

THE  SOLUTION  PROCESS  DESCRIBED  IN  GENERAL  TERMS 


The  main  task  is  the  evaluation  of  the  integrals  R.^ ,  Bg , 
and  Bg .  To  carry  out  the  discretization  needed  fur  this  purpose 
the  author  proposes  to  divide  the  planform  into  t"i angular 
subareas.  The  unknown  function  h  will  be  characterized  by  its 
values  at  the  corner  points  of  the  triangles.  For  the  interior 
of  the  triangles  an  interpolation  will  be  used.  The  corners  are 
numbered  by  two  indices  (i  and  j)  roughly  corren pan  i  i ng  to  the 
numbering  scheme  suggested  by  a  two-dimensional  coordinate 
system.  The  program  will  generate  for  each  of  these  corners  a 
table  which  gives  the  values  of  h  at  evenly  spiced  times.  It  is 
assumed  that  the  same  time  interval  At.  is  main*  line  1  throughout 
the  computation.  The  entries  within  these  tables  are  indexed  by 
a  subscript  k  (which  then  refers  fr  a  time  kAt).  In  the  memory 
the  function  h  then  appears  as  a  subscripted  variable  with  three 
subscripts,  the  first  two  for  the  corner  point  on  t ne  planform 
and  the  last  one  for  the  time. 

Assume  that  at  some  stage  in  l  ne  computation  Monte  tables 
have  been  generated  up  to  a  time  index  k .  The  above  equations 
then  allow  one  to  evaluate  ;t  h/  at  at  this  time  foi  >  1  1  stations. 

These  values  are  then  used  to  proceed  by  one  time  step.  The 

present  report  deals  with  the  eva]  nation  af  Jh/iu.  In  evaluating 
3h/9t  at  a  point  i,j  one  has  t<  e  mry  out  an  integral,  ion  over  the 
wing  and  the  wake  area.  This,  integration  appears  t  hen  as  a 
summation  over  the  points  of  the  wing  and  wake  ire a .  The 
summation  subscripts  will  be  m  irul  n .  .  The  listinction  between 

i  ,  j  on  one  hind  and  m  ,  r.  on  ‘Me  at. her  is  an  ii  oge.s  ,,  o  the 
distinction  between  the  eoor  1  i  naMm  x  ml  v  for  wii  h  Mie 
i  ntegro-di  f  f  ereiit  i  al  equal  i  or;  ,  s  t.  it  i  sf  i  .mi  Mo  umbral 

variables  of  i  nt.egr  it  i  on  mi  >,  whi  r,  o  ,  ;  ,  Ji;  •  .  •,  •  |-. 

integrals.)  The  boundary  'oni.ti  >i.  ,  ui-miv  :  m  ;;  v  is:;  ii  the 

point,  i  ,  j  eni,.-r  i  n  t.  be  ex;:r>- .  -n  m  ■  •  'nr  *  .  Ihe  point 

i,j  plays  a  si.mil  ir  '  ;  c.‘  •  .  : 


•M.-r 

-  j-  >. 


i  'j- 
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retardation  is  not  quite  constant  because  the  point  in  the 
field  can  be  regarded  as  a  substitute  for  a  whole  area  surround¬ 
ing  it . 

Assume  that  for  fixed  ij  and  fixed  mn  only  one  value  of  the 
retardation,  denoted  by  l  ,  is  obtained.  Then  Eq .  (7)  would 
still  be  applicable  but  all  coefficients  of  (mn)  l  would  be 

zero  except  for  l  =  I  .  However,  this  is  not  the  case  because  mn 
stands  for  an  area  in  the  £n  plane.  Moreover,  an  interpolation 
with  respect  to  time  is  needed  to  accommodate  nonintegral  values 
of  the  retardation  .  One,  therefore,  will  obtain  a  band  (in 
terms  of  8.)  of  the  values  C,.  ..  ,  ,  ,  in  which  these  quantities 

are  different  from  zero;  this  band  is  fairly  narrow. 

The  coefficients  C,.  .  „  arise  in  essence  by 

evaluating  the  integrals  3^ ,  Bg ,  or  Bg  at  fixed  ij,  separately 

for  the  individual  triangular  areas  or  for  the  hexagon 

surrounding  the  point  ij.  Adjacent  to  each  point  (mn)  there  are 

six  triangles.  Within  the  triangles  h  is  determined  by  its 

values  at  the  corners.  Accordingly,  the  integrands  for  one 

triangle  of  the  hexagon  depends  upon  the  values  of  h  at  the 

corner  points.  To  be  more  specific  one  evaluates  in  each  area 

the  integrands  at  certain  pivotal  points  and  multiplies  them  by  a 

factor  which  is  roughly  speaking  the  equivalent  of  the  area 

pertaining  to  the  pivotal  point.  The  contribution  of  one  pivotal 

(  2  )  (  3  ) 

point  depends  upon  the  value  of  h  and  h  and,  therefore, 

upon  the  values  of  h  at  the  corners  of  the  area.  One  obtains 

contributions  of  each  pivotal  point  to  the  values  of  C^. ^  (mn)  l 

for  those  values  of  mn  which  determine  h  wit.hin  the  area.  After 

the  procedure  has  been  carried  out.  for  all  ij  and  for  all  pivotal 

points  on  the  wing  and  the  wake,  the  evaluation  of  the 

coefficients  0.  .  is  complete.  The  report  discusses  in  detail: 
l  j  in  n 

1.  The  numbering  system  for  the  corners  of  the 
t  r i angles . 

?.  Tiie  integration  proeodur*-  in  general  . 


3.  The  representation  of  h  within  the  triangles 

and  within  the  areas  surrounding  a  control  point, 
and  transformations  required  to  obtain  well 
behaved  integrands. 

4.  A  criterion  to  decide  in  a  given  region  which 
precision  of  the  integration  formulae  is  needed. 
(Such  a  criterion  would  replace  a  decision  made 
otherwise  by  i nspections . ) 


V, 


SECTION  VI 


INDEXING  OF  CORNERS  AND  OF  AREAS 


The  numbering  scheme  is  obtained  in  the  following  manner. 
The  planform  is  first  divided  into  quadrangular  areas.  Triangles 
are  obtained  by  drawing  one  diagonal  within  each  quadrangle.  The 
corners  of  the  triangles  are,  of  course,  the  same  as  the  corners 
of  the  quadrangles.  The  corners  are  characterized  by  two 
indices,  the  first  one  for  the  (approximate)  chord  direction,  the 
second  for  the  span  direction.  The  subscripts  used  for  a  control 
point  are  i  and  j,  for  points  within  the  field  m  and  n.  The 
value  of  i  for  the  leading  edge  is  zero,  for  the  trailing  edge 


i  .  The  quadrangles  are  characterized  by  the  indices  at  their 


lower  left-hand  corner.  The  two  triangles  which  arise  after 
drawing  one  diagonal  within  the  quadrangle  are  distinguished  by 
further  subscripts  1  or  2  (Figure  3)*  It  may  appear  clumsy  to 
use  three  subscripts  to  identify  one  triangle  but  the  memory 
space  thus  taken  up  is  not  larger  than  the  memory  space  for  a 
numbering  by  one  index. 


The  geometry  of  the  problem  is  introduced  into  the 
computation  by  a  list  of  the  corner  coordinates  and  y^.  If 

the  wing  is  symmetric,  one  need  only  consider  the  right-half  of 
the  wing.  The  x  axis  is  then  the  line  of  symmetry.  The 
symmetric  and  antisymmetric  parts  of  the  motions  and  deformations 
of  the  wing  would  be  treated  separately.  We  write  down  the 
corner  numbers  pertaining  to  certain  areas. 


While  one  carries  out  computations  for  a  specific  triangle, 


the  indices  of  the  corners  are  denoted  by  and  n^.  One  has  for 


a  triangle  (mnl ) 


m^  =  m 


,  n  .j  =  n 


m^  =  m  +  1  ,  n,  =  n 


m^  =  m 


,  n  ?  =  n  +  1 
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(The  values  of  h,  E, ,  and  n  for  these  three  points  will  be  denoted 
by  h  ,  E  ^ ,  and  2.  =  1,  2,  3.  The  symbol  h  instead  of  h  is 

introduced  because  h  appears  already  as  subscripted  variable  hmnk 
in  the  program.  The  above  relation  ( Eq .  (8))  will  occur  only 
twice  in  the  treatment  of  a  specific  triangle,  the  first  time 
when  one  identifies  the  coordinates  of  the  corner  points,  for 
instance , 

m( 1 )  , n (  1  ) 

the  second  time  after  the  integration  for  the  area  has  been 

* 

carried  out,  when  one  identifies  h.  with  h_  _  .  Therefore 

_m(  2. )  ,  n  ( 2. ) 

there  is  no  need  to  introduce  the  mf  and  n^  in  the  program 
explicitly.  Furthermore,  one  has  for  a  triangle  (mn2) 

m1  =  m  ,  n1  =  n  +  1 

m2  =  m  +  1  ,  n^  =  n  (9) 

m^  =  m  +  1  ,  n^  =  n  +  1 

A  control  point  ij  appears  as  vertex  in  six  triangles 
surrounding  it.  For  the  whole  area  covered  by  these  six 
triangles  a  single  representation  for  h  is  used;  it  depends  upon 
the  values  of  h  at  the  corners  of  all  these  triangles,  including, 
of  course,  h^. 


# 

In  these  expressions  we  have  written  m ( l )  instead  of  m^  in 
order  to  avoid  two  level  subscripts.  In  general  we  shall  prefer 
notations  by  subscripts. 


While  one  deals  with  the  area  surrounding  a  specific  control 
point  ij,  these  points  are  numbered  from  1  to  7  (Figure  4).  One 
has 


m1  =  i 
m2  =  i 
m^  =  i  +  1 
m4  =  i  +  1 
m5  =  i 

S6  =  i  -  1 

*  i  -  1 


n2  =  j  -  1 

n3  =  j  -  1 

n  4  =  J 

n5  =  j  +  1 

n6  =  J  +  1 


(10) 


For  control  points  on  the  line  of  symmetry  (j  =  0)  only 
three  of  these  triangles  appear  in  computation.  Therefore,  the 
representation  within  the  area  depends  only  upon  five  values  of  h 
(see  F igure  5) . 


m1  =  i  ,  n1  =  0 

m^  =  i  +  1  ,  n2  =  0 

=  i  ,  n^  =  1  (11) 

m4  =  i  -  1  ,  n  4  =  1 

=  i  -  1  ,  n^  =  0 

The  function  h  appears  in  the  program  in  the  form  hmnk 
where  m  and  n  identify  the  corner  point  of  the  net  and  k  refers 
to  the  time  argument  (t  or  t).  Let  At  be  the  time  interval  for 
which  the  computation  is  carried  out.  Then 


t  =  kAt 


During  the  process  of  integrating  over  a  specific  subarea  the 
corner  points  are  characterized  by  one  subscript  (usually  l) .  We 
use  the  notation  h^.  If  the  i/At  is  not  an  integer,  we  then 


20 


we  denote 


write  h^(x).  If  we  want  to  use  the  argument  t  in  hmn . 
the  function  by  hmn,  so  that  for  x  =  kAt 

\m(t)  ’  ’  hmnk 

The  values  of  m  and  n  pertaining,  for  a  specific  area,  to  the 
index  2,  are  listed  in  Eqs.  (8),  (9),  (10),  and  (11).  Then 


VT) 


h 


mnk 


if 


is  an  integer. 


x/At  =  k 


SECTION  VII 


THE  INTEGRATION  PROCEDURE 


The  approximation  of  h  within  a  certain  area  will  appear  in 
the  form 


hU.S.n)  =  [4^  (  5  ,  n)  .  n)  .  .  .  ]A 


Here  the  known  functions  ,  4>2>  etc.,  are  the  elements  of  a  row 

vector.  A  is  a  matrix  which  depends  upon  the  geometry  of  the 

area  under  consideration;  h. ,  h2,’*'  the  elements  of  a  column 

vector,  are  the  values  of  h  at  the  corners  of  the  area.  The 

dimensions  of  A  are  naturally  determined  by  the  dimensions  of  the 

row  and  column  vectors.  In  most  cases  the  functions  <j>  ,<j>2...  are 

very  simple.  The  size  of  the  vectors  (which  depend  upon  the 

problem  in  question)  is  small.  One  has,  because  the  derivatives 
(2)  (  3 ) 

h'  and  h  are  formed  at  constant  t, 


These  expressions  are  now  substituted  into  the  integrals 
occurring  in  Eq .  (6),  and  the  integrations  are  carried  out  over 
the  individual  areas.  Taking,  for  simplicity,  only  part  of  the 
expression  Bg  one  then  has  to  evaluate 


d£d  n 


The  integration  is  to  be  carried  out  for  a  subarea  of  the  wing. 

In  most  cases  one  or  several  transformations  of  the  independent 
variables  are  necessary;  sometimes  the  functions  ^  are  expressed 
in  coordinates  different  from  the  system.  Further 
transformations  are  needed  in  order  to  generate  a  smooth 
integrand  in  an  area  in  which  the  independent  variables  range 
from  -1  to  +1.  (Triangles  in  the  interior  of  the  wing  at  a 
distance  from  control  point  ij  and  triangles  in  the  wake  are 

exceptions.  There  one  uses  the  form  B„.  The  representation  for 

( 2 )  ( T )  ' 

h  and  h  is  then  so  simple  that  one  can  integrate  directly 

in  the  plane,  using  special  formulae  for  triangular  areas.) 

To  fix  the  idea  for  the  other  cases  let  the  final  independent 

variables  be  p  and  q.  Then  one  obtains  for  the  above  integral 

+  1+1 

-1-1 

f (  3)  U,  n)  [<t>1  (  n.  <t>2(  n-  •  •  ]A 


The  variables  E,  and  n  are  expressed  in  terms  of  p  and  q. 
Therefore,  the  whole  integrand  appears  as  a  function  of  p  and  q. 
It  can  be  evaluated  by  Gaussian  integrations,  separately  for  the 
p  and  q  directions.  In  a  familiar  manner  one  then  chooses 
pivotal  points  and  multiplies  by  the  pertinent  weights  (say  w^ 
and  wq ) •  Let 


3  ( n) 
aCpTqT 


{f^2)  (E,,n)[4>i  ,  e 


1 

J 


♦  f  (  K  *  n )  [  4>  i  •  •  •  3 )  A 

•  i  n  y  n 


[  ip  1  ,  i>2 . . .  ] 


This  is  a  row  vector  of  the  same  dimension  as  the  column 


vectors 


h  i  (  t  ) 
h2  ( t ) 


Then  the  contribution  of  this  pivotal  point  to 


B9  is 


h^  (  x  )  +  ^ T  ^  + 

The  argument  x  is  given  by  t  -  ret ( £  -  x,  n  -  y )  where 

ret  =  (  1  -  M2)"1  [MU  -  x)  +  p] 

?  ?1/2 
p  =  [  (  5  -  x )  +  (  n  -  y )  '  J 


The  value  of  ret  must  be  evaluated  for  the  pivotal  point  under 
consideration.  In  the  above  expression  the  value  of  ret  is  the 
same  for  all  functions  h.,,  h2,  etc.  Let 

ret  =  ret/At 

and  ret+  and  ret  be  respectively  the  smallest  integer  greater  or 
equal  to  ret,  and  the  greatest  integer  smaller  than  ret.  Of 
course 


ret  +  -  ret  =  1 

Now  h^(x)  is  obtained  by  linear  interpolation  between 

h ^ ( k  -  ret+)  and  h^(k  -  ret  ) 

To  be  specific 

h  (x)  *  h?(k  -  ret+)(ret  -  ret  )  -  h ^ ( k  -  ret  )(ret  -  ret  ) 


Thus  one  obtains  as  contribution  of  one  pivotal  point 


4^  (ret  -  ret  -  ret  +  )  -  4^  (ret  -  retf)h1(k 

♦  ^(ret  -  ret  )h2(k  -  ret  +  )  -  ip  2  ( ret  -  ret+)h0(k 


ret  ) 
ret  ) 


Now 


To  obtain  the  form  of  Eq . 

4^  ( ret  - 


h.  = 

m  (  H  )  n  (  l ) 

(7)  one  accumulates 


ret  )  i  n  C 

i  j  .m^ ,n  ,retf 


and 


~  \p  ( ret  -  ret  *)  in  C 

i  j  ,mt ,n l , ret 


(This  is  the  second  occasion  where  the  expressions  mJ  are  used.) 
The  procedure  is  the  same  for  triangles  in  the  interior  of  the 
wing  and  not  adjacent  to  the  point  ij  and  triangles  of  the  wake, 
except  that  one  can  use  special  formulae  for  the  integration 
within  triangles. 


SECTION  VIII 


REPRESENTATION  OF  h  FOR  DISTANT  TRIANGLES  AND  TRANSFORMATIONS 
REQUIRED  FOR  THE  INTEGRATION  PROCESS 

The  determination  of  the  constants  C,.  ,  requires  do¬ 
ll  j  )  ,  (mn) , k  M 

loops  over  i,j,  m,  and  n.  It  is  probably  preferable  to  choose  m 

and  n  a3  variables  in  the  outer  do-loops,  because  this  reduces  the 

amount  of  temporary  storage.  For  triangles  not  having  the  control 

point  i,j  as  one  of  their  corners  and  not  adjacent  to  the  leading 

or  trailing  edge,  h  within  the  triangle  is  approximated  by  linear 

interpolation  between  the  values  at  the  corners.  The  coordinates 

of  all  corner  points  x  and  y  are  initially  stored  in  memory. 

mn  mn 

Then  the  program  proceeds  in  the  following  steps.  For  an  inner 
triangle  (mnl  or  mn2)  one  extracts  from  memory 

=  x(mU),  n(l)) 

8.  =  1,2,3  (13) 

«  y(mU),  n  U  ) ) 

The  values  m ( l )  and  n ( l )  are  found  in  Eqs.  (8)  and  (9).  It  is  not 

necessary  to  introduce  in  the  program  mil)  and  nil)  explicitly. 

We  write 


Let 


h_  _  (t)  =  h(x) 

m(i),n(t) 


l 
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The  approximation  for  the  triangle  in  question  is  written 


h  (  t  ,  f  ,  n ) 


Postulating  that  h  assumes  the  values  h1 ,  h0 ,  and  h^  at  the 
respective  corners  one  obtains 


In  principle,  at  least,  one  sets  up  the  matrix  M  and  then  forms 

- 1  (?)  (  a ) 

M  .  Actually  one  needs  only  h  c  =  3h/3f,  and  h  -  9h/3n.  One 

has 


Since  the  first  element  of  the  row  vector  on  the  right  is  zero  in 

both  cases  the  first  row  of  the  matrix  M  1  is  not  needed.  We 
write,  accordingly 


Because  of  the  simplicity  of  this  problem  one  can  write  down  the 
elements  of  A  directly  (without  calling  a  m  itrix  inversion 
routine).  On**  has 


D  e  t  -  L 1  fij  -  t  ^  rj 


(16) 


The  area  of  the  triangle  is 


Dot /P 


1 1 

-(A,2  * 

Au> 

1  2  * 

n  -,/Det 

i  3 

-n,/D e  t 

21  = 

-{A2P  + 

»;,3) 

22  ^ 

-  6^/Det 

A,,  =  l  ,/Del 


flow  one  applies  the  integration  formulae  for  *  r  i  angles  reproduced 
in  Appendix  D  of  Reference  1  and  determiner,  the  con*,  r  i  hut  ions  to 
C,  .  .  .  ,  .  _  according  to  the  procedure  in  Sect  ior,  VI. 


Triangle 

s  adjacent  to 

the  1c 

ad  i  ng 

-dge 

ir  (' 

'•  n , i  r  iet  e r  i  ?. 

od  by 

the  indie e s  ( m 

-  0  ,  n  ,  k  )  ,  k 

*  1  ,  . 

F  o  r 

t  h-sn 

t  lie 

repr  esent  at 

i  o  n  o 

h  are  found  in 

a  coordinate 

sys' em 

u  ,  v  i 

n  w  h  i 

ch  ‘ 

h"  u  direct 

ion  . 

normal  to  t h e 

leading  edge. 

and  f. 

9.  ’  ' 1  v 

am 

act  ■!  de  f  i  n 

•  ■d  by 

F  q  s .  (13)  and 

(  1  *0  .  To  t  r a n 

s  f  u  r’m 

to  the 

u  ,  v 

sys- 

"r  wt-  lid’ :  n 

(  F i gu re  6 ) 


'j 


KT 


-.■'V'vxv 


-  .  2  -2,1/2  -  .,-2  -2.1/2 

cos  a  -  n^/(£^  +  n^)  sin  a  =  ^/(^  +  n3' 


Here  a  is  the  local  sweep  angle.  Only  cos  a  and  sin  a,  not  a 
itself,  will  be  needed.  Now  one  sets 


t  =  u  cos  a  +  v  s  i  n  a 


n  =  -usina  +  v  cosa 


u  ^  cos  a  -  n  sin  a 


v  =  C  sin  a  +  n  cos  a 


3U,n)  .  . 
3  ( u  ,  vT 


At  the  leading  edge  one  has  h  =  0.  The  steady  state  problem 

1  /2  * 

suggests  that  h  -  u  .  One  therefore  has  for  a  triangle 
(  m  -  0  ,  n  ,  1  ) 


hne  has 


W2  ~1/2  -  .  , 
h=u  Up  hp  (  t  ) 


and  one  computes 


U1  "  u3  *  0 


u ,  »  cos  a  -  sin  a 

c.  r  t- 


The  P  r  -indt  1  -  GJ  nuert  distortion  changes  the  sweep  angle  and 
the  normal  to  the  leading  or  trailing  edges.  The  normal  in  the 
xy- system  Iocs  not  transform  into  the  normal  in  the  (distorted) 

f, r'- sys  t  ••m  .  Therefore,  we  have  ret  arned  in  the  present  report  to 
»  he  original  coordinates. 


-  -  -J \  ■  ‘  •  *  n.*  *-  *  *-  *-  ' 


V  V  ^  1 


•'  -T,jr-'r-' 11  '•' » -  ^r>g  v.1  VI  *rmj  t  f.'VTvj^:  w.-vj  m  Tvy.v.wy.  ».  „  ,..  ».„  ^  r.  n..  w  wwjw 


To  carry  out  the  integrations  needed  in  By  for  a  triangle  (m 

-  0,n,1)  one  has  to  form  h  and  h  J  .  One  obtains  from  Eq .  (19) 
and  (20) 


(2)  1  ,,-1/2 


-1/2  r 


(3) 


~2  u  cos  a  '  h 2 ^ T  ^ 


1  -1/2  .  „  -1/2  r  ,  « 

•^u  sin  a  Uj  h  2  (  t  ) 


(These  expressions  have  the  form  of  Eqs.  (12),  but  all  vectors 

(  2 )  (  3  ) 

have  dimension  1.)  The  functions  f  and  f  are  regular  within 

(2)  (3)  -1/2 

the  triangle,  but  h  and  h  for  u  =  0  behave  as  u 

Moreover,  one  must  take  the  triangular  shape  of  the  element  into 

account.  For  the  latter  purpose  we  introduce  a  system  of 

coordinates  centered  at  the  point  2  (Figure  7). 


u  -  -(u  -  u 2 ) 


v  =  -  (  v  -  v  2 ) 


u  «  U2  -  u  (rotation  by  tt  and 
translation) 


v  -  v2  -  v 


and  set 


v 

w  *  — 


Then 


3Un) 

3  ( u  ,  w  ) 


u  du  dw 


One  has  for  the  corner  points 


3  1 


AAii>  ^  j 


0 


u1  = 


(because  =  0,  and  U£  -  Moreover 


-  -<v,  -  »j) 


;3  -  -(*3  -  »2> 


v2  -  0 


The  limits  of  integration  are  then 


u  =  0  and  u  =  u? 


and 


w  =  w. 


<v3  ‘  V 


and 


w  =  w. 


(v,  -  Vj> 


To  eliminate  the  singularity  caused  by  the  factor 


u  1/2  -  ( u9  -  u)  1 /2 ,  we  set 


(u_  -  u)  =  pc 


-du  =  Tpdp 


3? 


The  limit  u  =  0  give  p  =  ;  the  limit  u  =  gives  p  =  0.  The 

Jacobian  becomes 

’  -2P“<1"<1P 

The  Jacobian  is  negative,  but  this  is  compensated  by  the  fact  that 
the  upper  limit  in  p  is  smaller  than  the  lower  limit. 

By  these  transformations  the  area  integration  becomes  a 
rectangle  in  the  wp-plane,  and  the  integral  is  a  smooth  function 
because  the  factor  p  in  the  Jacobian  cancels  the  denominator 

u1/2  »  ( u2  -  u)1/2  =  p.  This  integral  offers  no  difficulties  in 

carrying  out  a  Gaussian  integration  for  the  p  and  w  coordinates. 
One  may  introduce  a  further  transformation  so  that  the  region  of 
integration  extends  in  both  directions  from  -1  to  +1. 


Then 


w  -  (w1  +  ) /2  *  w(w1  -  w^)/2) 

P  =  (0  +  Up1 72  )/2  +  p(0  -  u^72  )/2 


3(5. n) 

3  ( p  ,  w  ) 


p  u  (  (  W1 


w3)/2)  -u^72 


One  chooses  the  pivotal  points  for  the  desired  order  of  Gaussian 

integration  in  the  p  and  w  integration,  and  determines  for  the 

pivotal  points  in  turn,  p,  w,  u,  v,  u,  v,  and  £,n.  Then  one 
evaluates  the  integrand  (including  the  Jacobian)  and  the 
retardation  for  the  pivotal  points  (and  in  the  case  of  Gaussian 
integration  multiplies  by  the  weights).  With  these  integrals  one 


determines  the  contributions  to  c(ij)  (mn)  ^  according  to  the 
procedure  of  Section  VII. 

For  the  triangle  (m  =  0,n,2)  one  introduces  again  the 
Quantities  defined  in  Eqs.  (13)  and  (14).  The  sweep  angle  a  is 
again  determined  by  the  above  formulae,  but  one  might  also 
introduce  auxiliary  points  4  and  5  (Figure  8)  and  introduce 


cos  a  =  ( n5  -  ni|)/[(^5  -  54)2  +  ( n5  ••  n4)2]1/2 
sin  a  =  U5  -  £;4)/[U5  -  £4 ) 2  +  (  n5  -  n4)2]172 


(21 


For  a  straight  leading  edge,  there  is,  of  course,  only  one  angle 
a.  Now  one  sets  for  a  triangle  (m  =  0,n,2). 


h(T ,£  f  n) 


Then  by  matching  h  at  the  points  2  and  3 


Then 


h(  ,  n,  t  ) 


[uc,v]A 


h^  (t ) 


Hence 


II 

00 

x: 

_  1  /p 

( 1 /2 ) u  cos  a,  v  sin  a 

A 

h?  (  t  ) 

_ 

h^  ( T  ) 

-  M 

(22) 


r  -i/2 

-(1/2 )u  sin  a,  v  cos  a 

A 

hp  (  x  ) 

L  J 

h  ^  (  t  ) 

To  take  the  triangular  shape  of  the  element  into  account  a 
new  system  of  coordinates  is  introduced,  so  that  the  side  2,3  of 

the  triangle  is  a  line  u  =  constant.  For  this  purpose  define 
cos  a  =  ( v  ^  -  v2)/[(u^  -  u2)2  +  (v^  -  v2)2]1/2 

sin  a  =  ( u  ^  -  u2>/[(u^  -  i^)2  +  (v^  -  Vp)2]172 

set 

u=ucosa+vsina  ,  u=ucosa-vsina 

v=-usina+v  cos  a  ,  v=usina+v  cos  a 


'-?1! 


and  compute 


I 

t 

t 

!  U2’V2  ’  U3,V3  ’  U2  =  U3  ’  and  v 2 

!  Furthermore,  introduce 


w  =  v/  u 


and 


compute  w,. 


and 


Then 


3U,  n) 
3 (u , w  ) 


u  d  u  dw 


Then 


(2) 


1  ~  —  1  / 2  -  -  -1/2 

7T  u  (cos  a  +  w  sin  a)  cos  a,  sin  a 


( t  ) 

Lh3(,) 


(3) 


[' 


1  -  - 1  /2  /  -  -,-1/2 
^  u  (cos  a  +  w  sin  a) 


1 

h^  ( t  ) 

ot 

A 

- 

h^  ( T  ) 

Actually,  one  simply  evaluates  u  and  v  for  the  pivotal  points  and 
substitutes  in  Eqs .  (22),  but  the  last  equation  shows  the 

singularity  that  arises  for  u  =  0. 

The  limits  of  integration  are  now 

u  =  0  and  u  =  u^ 

and  w  =  Wp  and  w  =  w^ 

-  -  i  /  p 

To  eliminate  the  singularity  caused  by  the  factor  u  -  we  set 


36 


=  it-1  ,  n  ^  =  n  +  i 


m4  -  it  •  n4  “  n  +  1 


m5  =  Xf1  ’  n5  =  n 


One  determines  ^  and  (  H  =  1...5) 


t’l  ’  ^2 


l  =  1  ...  5 


n£  ~  n2 


Furthermore,  one  evaluates  for  the  sweep  angle. 


oos  a  =  ^  /(ll  ♦  n^ ) 1 72  ,  sin  a  =  lh/(ll  ♦  n2)172 


-4/  v  ^4  T  "4. 


carries  out  the  rotation  of  the  coordinate  system 

u  =  f,  cos  a  -  n  sin  a  ;  £,  =  u  cos  a  +  v  sin  a 
v  =  £  sin  a  +  n  cos  a  ;  n  =  -u  sin  a  +  v  cos  a 

and  determines  the  values  of  u  and  v  for  points  1  through  5 
the  wake  there  is  no  singularity  of  h  at  the  trailing  edge. 


h  =  [1 


a- 

,u  , v  ]  b 

-C. 


,* 1 v  ’avv. -v  --  -■  /. 


The  expression  for 

a  term  of  the  f orm 
continuous  at  the 


To  determine 

at  the  points  2,  3 
expressed  by 


with 


Compute 


h  at  the  wing  differs  from  that  for  the  wake  by 
3/2 

(-u)  .  Then  h  and  its  first  derivatives  are 

trailing  edge.  Therefore  on  the  wing 

h  =  [1,u,v,(-u)3/2] 


a" 

b 

c 

d 


h  in  the  triangle  (m  -  i  -  1  ,  n  ,  2 )  ,  it 


A  - 


i  n 

add  i  t 

ion 

h2 

*a 

=  M 

b 

h5 

c 

_h3_ 

d 

0 

0 

0 

10  v , 


1  U5  V5 


J  U3  V3 


0 

0 


( -U^) 


3/2 


s  matched 
i  s 


- 1 


A  =  M 


h 


[ 1  , u , v , 0]  A 


1 1 
h 
h 


4 

5 
3 


in  the  wake, 

and 

h  =  [1,u,v,(-u^/2]A 

on  the  wing.  The  representation  for  the  wake  can,  of  course,  be 
obtained  more  directly,  but  the  present  form  is  somewhat  easier  to 
program . 


si 


For  the  triangle  (m  =  i  ,n,1)  one  matches  at  the  points  1, 
2,  3,  and  5.  Then 


h  =  [ 1  ,  u , v  ,  ( -u ) ^ /? ] 


hr 


and  the  pertinent  matrix  M  is  given  by 


C  om  pate 


A  =  M 

Then 


h  =  [1 ,u,v,(-u}3/2]A 


Now  one  has  for  the  triangle  (m  =  ,n,l ) 


h(2)  =  jcos  a[  0  ,  1  ,  0  ,  -  |(-u)1/2]  *  sin  «[0,0,1,0]\A 


h^  =  |  sin  a[  0  ,  -  !  ,  0  ,  |(-a)1/<‘]  +  cos  ot[  0,0,  1  ,  0  J  ^  A 


and  for  the  triangle  (m  =  i 


,  .  n  , .?  ) 


h(2)  =  Jcos  a[0,  1  ,0,-(3/2)  (-u)  1 /2]  ♦  sin  a[  0 , 0  ,  1  ,  0  ]  j  A. 

(  7)  (  1/2  / 

hv:i;  =  |  si  n  a[0, -1  ,0,  (  3/2)  (-u)  ]  +  cos  a[0,0,l,0jjA 

(Of  course  with  different  matrices  A.) 

The  further  steps  of  the  treatment  are  different  for  the 
triangles  (m  =  i t - 1 , n , 1 )  and  (m  =  it~1,n,2).  The  side  (1-3)  of 

the  triangle  (m  =  i  t - 1 , n ,  1 )  is  not  necessarily  parallel  to  the 

trailing  edge,  therefore,  a  rotation  of  the  coordinate  system  is 

carried  out  so  that  u1  equals  u^  (Figure  10). 

cos  a  =  ( v  .j  -  v3)/[(u1  -  u^)^  +  ( v  1  -  \'3)2]1/2 

sin  a  =  (u1  -  u3)/[(u1  -  u^)^  +  (v1  -  v3)2]1/2 


u  =  u  cos  a  +  v  sin  a  ,  u  =  u  cos  a  -  v  sin  a 


c 

o 


et 


v=-usina+vcosa  ,  v=usina+veosa 


w  =  v/  u 


Then 


The  limits  of  integration  are  0  and  u^  for  a  and  w^  and 

1  /? 

for  w.  The  factor  (-u)  gives  rise  to  a  term 

,  -  A  /2  ,  s  1 /2 

(u)  (-cos  a  +w  sin  a) 

Therefore  we  introduce 

__  O 

u  =  p4- 


Then 


9  (  E, ,  n ) 
3  fp  ,  w  5 


2up 


Again  we  transform  p  and  w  so  that  in  either  case  the  1 
are  -  1  and  + 1  . 


p  =  (  (u  1  /2  )/2)  -  p(  (u  ^  '2  )/2) 


w  =  ( w  *  w1  )/2  ♦  w  ( w  1  -  w?)/2 


Then 


=  (1/2)up(u  )1/2(w  *  w  ) 

3  (  p  ,  w  )  3 


The  argument  of  h^_  needs  special  attention.  One  has 


hr  (  t  ) 


h ,,  (  x 


—  (  r 


C2)  ) 


with 


4  d 


imi  ts 


T 


ret 


Therefore,  one  forma 


ret  -  [ret  +  rt.  -  f„)]/At 


deter  mi  nes  "he  values  of  f.  from  ret  and  ret  ,  and  accumulates 

+  — 

f  hose  cent  r  i  but  i  ons  in  h.,  for  t  -  ret  md  ret  ,  Actually,  thi 
i  the  prone  lure  for  all  points  in  the  wake. 

Tor  the  triangle  (m  =  i,  ,  n ,  ? )  one  first  shifts  the  eoor 
1 1  ri  at  •?  s  v  .>  t  em  to  poi  nt  1  (  F  i  gure  11). 


u^  ,  u  -  u  ♦  u 


v  ■  v  -  v^  ,  v  -  v  ♦  v 


W  =  V/  u 


_  J  /2  .  ,  -  ,  1  /. 

The  t  erm  ,  -  u  becomes  (  -  u  ^  ♦  u  ) 


T o  count er  ict  t ht 


■iinviiari)  y  caused  by  t  his  term  we  set 


So  far  one  has  the  integration  over  a  rectangle.  To  rehuee  this 
to  the  integration  over  a  square  we  set 


,  ,  J/2  J  /2 

p  =  (  (  _  u  ^  )  /2  )  -  p(-u^)  /  2 


w  =  (wQ  +  w^)/2  +  w ( w ^  -  w2 ) / ? 


3  (  Sh )  /  «  /t  n  -  /  1 /2  ,  , 

-  _  — =  (1/2)pu(-u  )  (wh  -  w0) 

3 (p , w  )  ^ 


SECTION  IX 


APPROXIMATION  FOR  h  IN  THE  AREA  SURROUNDING  A  CONTROL  POINT 
AND  TRANSFORMATIONS  REQUIRED  FOR  THE  INTEGRATION 

The  right-hand  side  of  the  integral  equation  (that  is  the 
expression  Bg)  is  particularly  sensitive  to  the  approximation 
chosen  for  the  function  h.  In  the  computation  only  values  of  h  at 
the  corner  points  of  the  grid  are  available.  A  local  mesh  refine¬ 
ment  cannot  be  carried  out,  because  eventually  each  corner  point 
will  be  a  control  point;  ultimately  one  would  have  an  overall 
refinement  of  the  grid. 

The  fact  that  h(2)(x.,y.)  =  f£(x.,y.)  and  h'-;(x  ,y  )  = 

11  d  X  1  l  11 

yy(x.,y.)  occur  in  a  number  of  expressions  contributing  to 
3h(x^,y.)/dt  may  prove  a  source  of  uncertainty.  The  contributions 
of  some  of  these  terms  do  not  go  to  zero  as  one  considers  a 
smaller  and  smaller  neighborhood  of  the  control  point,  (which  shows 
that  these  terms  are  indeed  important).  There  may  be  a  partial 
cancellation  of  these  terms.  Such  cancellations  would  result  in  a 
loss  or  accuracy,  if  one  approximates  in  different  triangles 
adjacent  to  the  control  point  the  functions  h^*’^  and  h^'  by 
different  expressions.  It  is,  therefore,  preferable  to  use  one 
analytic  expression  for  h  throughout  the  entire  area  surrounding 
the  control  point  under  consideration.  Actually,  one  cannot 
guarantee  that  h  is  an  analytic  function  throughout  the  area,  and 
one  cannot  claim  that  such  a  representation  has  higher  accuracy 
that  an  approximation  by  piecewise  linear  function.”'.  The  higher 
order  approximation  is  solely  used  to  avoid  an  inconsistency.  In 
the  realization  of  the  method  which  we  have  in  mind,  the  function 
w  (x ,y )  has  the  form  of  some  function  in  x  and  y  mul*i plied  by  a 
step  function  in  time  (details  in  Section  XT).  The  pv-ntial  then 
approaches  a  steady  state.  The  steady  state  is  given  by  a  smooth 
function  provided ,  of  course,  that  t  he  dcp-ni  i»n.v  f  w ,  upon  x  and 
y  is  smooth .  The  approx  i  mat  ion  of  h  in  t  hv  v  i  v  ■  ”  >.  mnt  rol 

point  by  a  polynomial  is  t  lion  t  i  vast  lg.-ius  b“c  »...••••  of  its  greater 
ae  uiracy  . 


In  studying  the  question  of  accuracy  in  detail,  one  can 
study  a  fictitious  problem  where,  in  the  vicinity  of  a  control 
point,  h  is  exactly  given  by  second  degree  polynomial.  (One 
evaluates  the  exact  contribution  to  the  upwash,  once  using  the 
exact  formula  for  h,  a  second  time  with  a  linear  approximation  for 
each  of  the  six  triangles  and  compares  the  results.) 

The  choice  of  the  approximation  for  h  in  the  area 
surrounding  a  control  point  requires  some  discussion.  A 
polynomial  of  the  second  degree  contains  six  arbitrary  constants 
but  one  must  match  h  at  seven  points  (six  corner:;  of  the  hexagon 
and  the  center).  A  third  degree  polynomial,  on  the  other  hand, 
contains  too  many  arbitrary  constants. 

The  situation  can  be  discussed  in  detail  for  a  regular 
hexagon,  and  nearly  all  the  results  can  be  carried  over  to  a 
•hexagon  which  arises  by  an  affine  t r ansf orma t i on . 

To  see  the  question  in  a  more  general  light  we  carry  out  the 
initial  discussion  for  a  regular  octagon.  Assume  initially  that 
the  corners  lie  on  a  circle  with  radius  1,  that  the  center  lies  at 
the  origin  of  the  x,y-system,  and  that  two  corners  lie  on  the  x- 
axis.  First  we  disregard  the  matching  at  the  origin.  Let 
arctg(y/x)  =  0.  The  matching  is  then  carried  out  at  equidistant 
values  of  6,  namely  0  =  ,  n  =  1...8.  It  can  be  accomplished 

by  a  discrete  Fourier  analysis.  The  functions  involved  are  1, 
cos  0,  sin  0,  cos  20,  sin  20,  cos  30,  sin  30,  and  cos  He.  The 
function  sin  He  does  not  appear  because  it  is  zero  at  all  matching 
points. 

The  case  in  which  no  corners  lie  on  the  x-axis  can  be 
studied  by  rotating  the  coordinate  system  by,  say,  an  angle  c. 

The  term  cos  He  then  appears  as 

e  o  s  '  H  (  0  -  c  )  )  =  c  o  s  (  H  t  )  c  o  s  (  H  0  )  *  s  i  n  (  H )  s  i  n  (  H  0  ) 

and  one  must  admit  sin  (Ho)  as  we l  1  as  cos  (Ho''.  Th  ■  excluded 


function  is  then 


sin(4(0  -  e ) )  =  cos ( 4c )si n( 46)  -  si n ( 4c ) cos ( 4 9) 

A  formulation  in  which  c  need  not  be  specified  is  obtained  if  one 
admits  cos  (40)  as  well  as  sin  (40)  and  postulates  that  the  sum  of 
the  square  of  their  coefficients  be  minimized. 

If  one  omits  in  this  problem  cos(4e)  and  sin(4e)  and 
postulates  that  the  sum  of  the  square  of  the  errors  be  minimized, 
then  it  follows  from  the  properties  of  a  finite  Fourier 
decomposition  that  all  coefficients  up  to  those  of  cos(30)  and 
s i n ( 4  0 )  remain  unchanged.  The  errors  are  those  caused  by  the 
omitted  term  cos (40)  or  cos (4(9  -  c)).  They  have  the  form  ±1 
multiplied  by  a  constant,  where  adjacent  points  have  opposite 
signs.  This  is  a  reasonable  approximation. 

Actually,  one  carries  out  the  approximation  in  terms  of  a 

polynomial  in  x  and  y;  in  the  case  under  consideration  this  is  a 

fourth  degree  polynomial.  The  constant  in  the  polynomial  is  the 

only  term  which  affects  the  value  at  the  center  of  the  octagon 

(the  origin  of  the  x , y- system) .  Postulating  a  perfect  match  at 

the  center,  one  has  to  match  the  vni  ies  of  h  -  h„.„,.  _„at  the 

corner  center 

corners.  For  this  purpose  one  has  the  following  expressions  at 
one's  disposal 

.?  ?  3  P  3  4  3  p  3  4 

x,  y,  x  ,  xy,  y  ,  x,  x  y,  xv  ,  v  ,  x  ,  xy,  x  y,  xy  ,  and  y 

The  expressions  eos(nO)  and  r.in(nO)  of  the  matching  by  a  Fourier 
decomposition  can  also  be  written  (along  the  circle  of  radius  one) 
as  polynomials  in  x  and  y.  They  are 

a  j  :>  ->  ?  q  a  e  3 

1  ,  x,  y,  x  -  y  ‘  ,  2xy ,  x~  -  3  x  y '  ,  3x‘  y  -  y  ,  x  -  bx^y"  +  yJ 

and,  is  a  counterpart  to  sin  4  0,  4x\y  -  4xy 

The  following  observations  are;  port!  nent  f  or  a  comparison  of 


error  the  contribution  of  the  omitted  third  degree  terms;  they  are 

of  the  form  ±1  multiplied  by  a  constant.  For  this  omitted  term 

3 

the  deviation  of  h  in  the  interior  of  the  hexagon  behave  as  r 
2  o  p 

(with  r"  =  +  y  ).  Accordingly,  deviations  are  rather  small  at 

some  distance  from  the  outer  edges.  To  allow  for  matching  at  the 
center,  one  also  admits  the  constant  1. 

We  ;..ve  assumed  in  the  derivation  that  one  has  a  perfect 

match  t  ,e  center.  For  the  regular  hexagon  this  does  not  affect 

the  minimization  process;  if  h  -  1  in  the  middle  an  1  zero  at  all 

o  2 

corner  points,  then  the  expression  1  -  x  -  y  gives  j  perfect 
matching.  The  requirements  for  perfect,  mat  oh  i  ng  by  tnir-d  degree 
polynomials  have  been  stated  a  bo  v*_* .  :  •  v-,n  b--  r;.;  r  that  for 

hexagons  that  do  not  deviate  :  oc  sv---np'  ,  from  mg.,  1  nr  hexagons 
the  conditions  developed  her’.-  a-.-  .-.v  ;  .-•••'  i  •  ry  . 

The  numbering  of  the  p  <  h-  x  •.  ••  n.  iinl  ,  ling 

the  control  point  itself  g  :■  •  ■  :  ‘  r  :■  ’  ’  .•  ••  •  4  .  V  :r 

every  control  point  one  i  :  *-•••..  ■  ■  •  ro¬ 

se  v  e  n  points  who  r e  h  i s  g i  .  -  -  . 


where  m  ( l )  and  n  ( l )  .ire  sef  .  •  ;  .  '  .  :••••  \  ■ 

point  in  the  i  nteri  or .  W  o  inf--  .  .  >  •  •••,  or 

at  the  cent  ml  po  1  r;t 


If  one  exp  res  s  03  h  by  a  qua  iraf  i  .-  fun  u  ion 
coefficients  at  one’.-  disposal.  I'M--  ippm 


expression  is 


:  u  l.  he  form 


t  fieri 

:'i  -i  <_  ; 


wr i tt cn 


2  2 

x  +  y  =  1 :  In  the  first  list  the  term  1  does  not  appear, 

2  2 

therefore,  the  terms  x',  xy,  y  of  the  first  list  express,  as  far 

as  the  matching  is  concerned  the  same  family  of  functions  as  1, 

x^  -  y'  and  2xy.  In  the  first  list  some  of  the  expressions  of  the 

third  and  fourth  order  can  be  expressed  by  terms  already  included 

2  2 

in  the  list  because  x"  +  y  =  1.  They  are 

2  33  222243  3 

xy  =  x  -  x  ,  y=y-yx,xy  =  x  -  x  ,  xy  -  -xy  -  x  y 


Accordingly,  they  ought  to  be  omitted.  Besides  one  has  the 

3  3 

expression  sin  40  =  4x  y  -  4xy  which  vanishes.  Accordingly,  the 

first  list  contains  a  number  of  terms  which  are  linearly 

dependent.  The  matching  by  a  fourth  order  polynomial  in  x  and  y 

is  therefore  best  carried  out  with  terms  of  the  second  list,  where 

2  2 

the  first  term  1  is  expressed  as  x~  +  y'.  Accordingly,  one  has 


and 


2  2  3 

x,  y,  x  ,  xy,  y  ,  xJ 

4x^y  -  4xy2 


3xy‘ 


9 

3x '  y 


r  2  2 

6x  y 


4 

+  y 


with  the  postulate  that  the  sum  of  the  squares  of  the  coefficients 
of  the  last  two  terms  be  minimized  for  the  regular  octagon.  This 
may  appear  nearly  trivial  for  the  regular  octagon  but,  by  using 
the  polynominal  expressions,  the  procedure  can  be  carried  over  to 
general  octagons.  Of  course,  the  nature  of  the  approximation  is 
then  only  loosely  defined. 


The  case  of  an  octagon  has  been  chosen  because  it  shows 

clearly  how  the  fact  that  x'  +  y1-  =  1  at  the  corners, 

eliminates  some  of  the  expressions  of  a  general  pol ynomi al  (because 

they  are  linearly  dependent  upon  expressions  already  included).  In 

the  case  of  a  hexagon  a  polynomial  of  the  third  degree  is 

sufficient.  One  then  admits  the  terms  x,  y,  x‘",  xy,  y  ',  x  - 
2  9  7 

3xy  ,  and  3x  -  y  (after  one  has  matched  at  the  center),  and 
minimizes  the  sum  of  the  squares  of  the  coefficients  of  the  last 
two  expressions.  If  one  omits  the  terms  of  the  third  degree  and 


mi  n  i  :n i  z e s 


as 


the 


urn  of  the  square;;  or  the  errors, 


one  obtains 


hU.n.k)  =  ( 1 


, S.n.C  ,  Sn.n  )  •  _ 

a^  (k  ) 


(Here  £  and  n  on  the  right  are  considered  as  functions  of  r  and  n 
(Eq.  23).) 

Matching  at  points  1  to  7  would  lead  to  the  requirement  that 


M  be  equal  to 


h  1  (  t  ) 


lVt)J 


where 


g  h  p 

r,3n3 


f’6n6 
f- 1 "  I 


Written  in  detail  this  would  require  that 


l  Muai  -  hi 
i  =  1 


i  *  1  ...  7 


be  zero.  In  general  this  is  impossible.  In  accordance  with  the 
discussion  at  the  beginning  of  this  section  we  impose  the  weaker 
r  equ  i  r ern  en  t 


Hence,  by  differentiation  with  respect  to  a 


Mt.  J,(Mua*  -  hi>  ■  0  =  ’•••« 


This  is  written  in  terms  of  matrices 


mm:  -  m  : 


H  en  ce 


with 


r  +  -  -  1  * 

A  -  [M  Mj  M 


This  procedure  is  known  as  forming  the  generalized  inverse  of  the 
matrix  M.  In  the  computation  one  first  builds  up  the  matrix  M  and 
subsequently  forms  M*M,  [M+M]  1 ,  and 

A  =  [MfM]_1M+ 


The  representation  for  h  is  then  given  by 


h  ( f, ,  n ,  k  )  =  [  1  ,  f  ,  n ,  ,  K  n ,  j  A 


Perfect  matching  at  all  seven  points  is  achieved  by 
admitting  expressions  of  the  t  hir'd  degree,  but  according  to  the 
above  discussion  the  third  degree  terms  are  chosen  in  a  special 

f  orm . 


As  we  postulate  perfect  matching  at  all  points,  we  introduce 


h  ( £, ,  n)  =  h  ( 5 ,  n)  -  h(  e,1  ,  n1 )  , 


that  is,  we  anticipate  the  matching  at  the  pivotal  point  1.  Then 
we  have  a  representation 


hU.n)  =  [£,n,C2,*;n,n2,U3  -  3^n2),(3C2n-  n3) 


Notice  the  special  form  of  the  terms  of  third  degree.  We 
partition  the  vector  Oa^-*] 

[a  1  .  .  . ]  =  [u1...u5  :  v 1  ,v2] 

This  means  that  we  denote  the  coefficients  of  the  linear  and 
quadratic  terms  by  u^  .  .  ,ur  and  the  coefficients  of  the  two  third 
degree  terms  by  v^  and  v^.  Then 


=  hU,n) 


We  choose  arbitrarily  five  corners  of  the  hexagon.  For  instance, 
the  corners  2  to  6.  Next,  define  a  5  by  5  matrix  B  by 


- 

~ 

U,nV2  Vn.n2] 

V 

♦  U3  -  3Cn2, 3C2n  -  n3  ] 

V1 

U5 

V  2 

B , 
B 
B 
B  . 


i  1 

“  V  +  1 

i  2 

=  Vi 

_  o 

i  3 

^i*i 

i  4 

Vi 

i  5 

“  nH 

i  =  1  ...  5 


n . 


a  5  by  2  matrix  C  by 


C i 1  "  ^  i  +  1  ^  3^i  +  1  n i  +  1 


C i  2  3Ci*1  ni  +  1  ”iH 


i  =  1  ...  5 


a  row  vector  o  with  5  elements  considered  as  a  1  by  ‘5  matrix  by 

-2  -  -  -2 
-  “  6  y  >  G  2  ~  hy>  C  ^  =  C  y  »  ^  J_)  =  o  y  n  y  »  f‘  =  h  y 

and  another  row  vector  d  with  two  elements 


--2  --2  -3 
d1  =  E,y  -  36ydy  and  d2  =  36yhy  -  hy 


Finally,  we  introduce  the  column  vectors  h  with  6  elements,  and 
h,  which  is  the  vector  h  truncated  after  the  fifth  element. 


h  . 
l 

=  h  ( 6  i  +  t  ,  n  i  + 1  ) 

,  i  =  1  ...  6 

h  . 
l 

h(6i+I ,oi+1 ) 

,  i  =  1  ...  5 

The  matching  conditions  at  the  points  2  through  6  and  at  point  7 
are  now  written  separately. 

B  u  +  C  v  =  h  (26 


■  u  +  d  v  = 


One  obtains  from  Eq .  (26) 

u  =  B”  (h  -  Cv)  (21 

Substituting  this  into  Eq .  (27)  and  collecting  terms  with  v  one 
obtains 


■V  ♦  fh  =  0 


with 


(30) 


e  =  d  -  cB  1 C 

f  -  [f, _ ,f6]  =  [cb-1  :  -  i ]  (3i ) 

(e  and  f  are  row  vectors  with  2  and  6  elements,  respectively.) 

Eq .  (29)  is  one  equation  for  the  two  unknowns  v]  and  v,;.  Written 
explicitly  it  reads 


6 

e.  v  +  e_v0  +  £  f-h.  =  0  (32) 

II  c.  t—  j  _  i  J  J 

Eq.  (32)  is  obtained  by  postulating  that  the  function  h(E,,n)  is 
matched  at  points  2  to  7.  As  Eq.  (32)  is  one  equation  with  two 
unknowns,  the  result  is  not  uniquely  determined.  We  now  introduce 
the  postulate  that  the  sum  of  the  squares  of  the  coefficients  of 
the  third  degree  terms  be  a  minimum. 

( 1 /2) ( v1 2  +  v?2 )  =  mi n 


Eq .  (32)  is  a  constraint  imposed  on  this  minimization.  Alter 
introducing  a  Lagrange  multiplier  \,  one  must  determine  the  values 
of  v 1  and  v2  for  which  the  expression 


( 1 /2  )  (  vi  2  -  v 2 2  ) 


6 

He,v,  .  .V2  •  fj  hjj 


is  stationary.  By  differentiation  with  respect  to  v1  and  v2  one 
obtai  ns 


v.  +Ae.  =0  ,  i=1,2 

l  l 


and  -iy  substitution  into  the  constraint  Eq.  (32) 


Hence 


,  2  2,-1  +  ~ 
v  =  -  ( e1  +  e2  )  e  f  h 

(One  remembers  that  v  is  a  column  vector  and  e  a  row  vector.) 
This  is  rewritten  in  the  form 


v  =  M^h 


with 

M2  =  -(e^  2  +  e22)  1  e  +  f> 


Then  from  Eq .  (28) 


u  =  M  ^  h 


with 


M,  =  B 


-1 


-  CM, 


Here  Ic  is  the  five  by  five  identity  matrix.  Then 
b 


a 


a 


=  M^h 


with 


In  the  program  one  generates  the  matrices  B  and  C  and  the  vector 
c  and  d  and  computes 


CD  L 


,  2  2,-1  V 

-(e1  +  e2  )  e  f 


8  p5  ;  °  -  cm2 


Then  one  has 


h(l,~r\)  -  U,n,£2,Cn,n2,  U3  -  SSn^M^n  -  3nJ)]A|: 


=2*  ,?2-  _  0-3' 


The  integration  procedure  is  now  described  for  the  approximation 
to  h  discussed  first,  (Eq.  25).  One  has  for  a  control  point  in 


the  interior 


[0, 1 ,0,2t,n,0]A  : 


h1  (t  ) 


h  y  (  t  ) 


[0,0,1 ,0,£,2n]A  : 


h  1  (  t  ) 


h  ^  ( i  ) 


The  expressions  to  be  evaluated  for  the  area  adjacent  to  the 
control  point  ij  are  Bg  and  Bg.  Bg  is  an  area  integral.  The 
integrals  are  evaluated  separately  for  the  six  triangles 
containing  the  point  (ij). 


Consider  a  triangle  with  corners  1,  H ,  l + 1 .  For  1=1 
identify  the  point  8.+1  with  the  point  2.  We  introduce  a  sys 
(u,v)  in  which  the  side  (  e. ,  £  +  1)  of  the  triangle  is  a  line 

u  =  const.  For  this  purpose  we  set 


cos  a  = 

<Vi 

^ 2.  ^  +  ^  nSl  +  1 

si  n  a  = 

(  .1 

( n£  +  1 

vP^I 

II 

C 

cos  a 

+  V 

s  i  n  a  , 

,  u  =  E,  cos  a 

2,1  /L 


n  -  -u  sin  a  +  v  cos  a  ,  v  =  F  sin  a  -  n  cos  a 

Furthermore . 


Then 


w  =  v/u  ,  v  =  uw 
3U.n)  .  ,, 


One  now  has  an  integration  over  a  rectangular  area 


0  <  U  <  Ut  ,  Wt  <  v  wt4l 


To  obtain  in  both  directions  limits  -1  and  +1  we  set 


u  =  u,/2  +  wUj/2  ,  w  »  (wl  +  1  +  wt)/2  +  w(wd  +  1  "  wi^/2 


Then 


3.(_§_L_nJ_  =  !  u  u  ( w  -  w  ) 

4  u  ’V'Vl  V 

3 (u , w  ) 


The  integrand  of  contains  h  in  the  forms 


h  ^  2  ^ ( i , f;  , n )  -  h  ^  2  ^  ( t  ,  x  ,  y  )  and  h  ^ ( t  ,  £ , n )  -  h  ^ '  ( t ,  x  ,  y  ) 


'*  n"'  .  .*  ■  •  /■  .  ,  ’ -\-v  . 


b3 


The  contributions  with  arguments  (t,x,y)  are  evaluated  separately 

(  2 ) 

but  for  the  same  pivotal  points  as  h  (t ,C,n) .  One  has 


h  ^  ^  ( t  ,x  .  .  ,y  .  .) 

1  j  1  j 


3)  (t  x  y  ) 
n  it,xij,yij; 


=  [ 0, 1 , 0, 0, 0,0]A 


=  [0,0,1 ,  0 , 0 , 0  ]  A 


h  1  ( t )' 


h^  ( t ) 


h1  (t ) 


h  ^  ( t ) 


For  these  expressions  the  retardation  is  zero.  These  contribu¬ 
tions  are  accumulated  in  C  _  _  ,  l  =  1...7. 

ij.mU)  ,  n(  l )  ,0 

The  expression  Bg  ( Eq .  (6))  is  a  contour  integral  around  the 
area  formed  by  the  six  triangles  adjacent  to  the  control  point 
(ij).  It  is  evaluated  separately  for  the  sides  ( l ,  S.  +  1)  of  the 
individual  triangles.  Because  point  1  is  given  one  has 


’  X i  :  ’  n  '  yij 


n-y-n  ,  ?  -  x  -  C 


Again  the  transformation  ( Eq .  (33))  is  carried  out.  The  side 
(l,  SL  +  1)  is  a  line  u  =  const.  Therefore, 

d£  «  d£  =  dv  sin  a  ,  dn  =  n  »  dv  cos  a 


,  2  2, 1/2 

p  =  (u^  +  v  ) 


The  variable  of  integration  in  Bn  is  now  dv,  the  limits  are  v, 


and  v.  . .  By  setting 


wl  =  ^1  +  1  +  v£,^/2  +  v  (  v  ^  +  -j  -  v  ^  )  /  P 


one  obtains  -1  and  +1  for  the  limits  in  v,  and 


JTtS 


V  vm  vp  1  w"  V-  V*? rjtvrj  ww  *  i  »  v*  v  JVWVV  ’/v  ."TV*  r*  '  *  ■'■» ' "  ’  v  \  ,v  ', 


* 

i 

/ 

*• 

-* 


w 

I 

/ 


C 

i 


h 


► 

h 


d  v 
dv 


( v 


i  + 1 


V/2 


For  an  area  adjacent  to  a  control  point  next  to  the  leading 

edge,  one  proceeds  as  follows.  At  the  leading  edge  the  potential 

is  zero;  there  no  control  points  are  needed.  The  solution  for 

1  /  2 

steady  flow  suggests  that  the  potential  has  a  behavior  as  u 
where  u  is  a  coordinate  normal  to  the  leading  edge.  (This  idea 
had  already  been  used  for  distant  triangles).  For  control  points 
with  i  =  1  (next  to  the  leading  edge)  we  then  proceed  as  follows. 
One  determines 


£n=x_  _  ,  n0=y_  _  £  =  1  ...  7 

m(0,n(0  m  ( Jt )  ,  n  ( {, ) 

with  and  n  by  Eq .  (10)  (Figure  i|)  .  Now  the  origin  is  shifted 
to  point  7. 


l 


n 


l 


=  n, 


<L  =  1  ...  6 


Next,  a  Cartesian  uv-system  is  introduced  with  the  u-axis  pointing 
downstream  and  normal  to  the  leading  edge. 

-  2  -2.1/2  -  .  2  -  2.1/2 
cos  a  =  n6/U6  +  hg  )  ,  sin  a  =  ^/(£.ft  ♦  n6  ) 

E  =  u  cos  a  +  v  sin  a  ,  u  =  C  cos  a  -  n  sin  a 

n  =  -usina  +  vcosot  ,  v  *  5  sin  a  *  n  cos  a 

One  determines  u  and  v  for  the  seven  points  at  which  one 
wants  to  match  h^  at  those  two  points  which  lie  at  t. he  leading 
edge  h  *  0 .  his  now  represented  in  the  form 


h  <  r  ,  n )  *  u  1  [  1  ,  u  ,  v  ,  a  v  ,  v  ‘  | 
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In  the  row  vector  an  element  u  has  been  omitted  for  the  following 
reason.  For  a  configuration  in  which  the  lines  2,5  and  3,^  ate 
parallel  to  6,7,  one  can  find  an  expression 

b1  +  b2u  +  b^u2 


which  vanishes  for  u  =  u,  =  u„  =  ur  and  u  =  u_,  =  u,,  with  coeffi- 

12  5  3  a 

cients  b^ ,  b^.  and  b^  different  from  zero.  If  one  would  include 

the  term  u2  ,  one  would  obtain  a  degenerate  matrix.  If  the 
configuration  is  close  to  the  one  just  described,  the  matrix  would 
be  close  to  degeneracy  with  the  consequence  that  some  terms  in  the 
inverse  might  be  very  large. 

In  this  case  perfect  matching  is  possible.  The  matching 
conditions  are  expressed  as  follows.  One  postulates 


with 


f-  1 

r  -i/2-| 

a , 

h ,  u , 

1 

1  1 

. 

- 

• 

• 

*  -1/2 

a  r 

h,  u 

5_ 

5  6 

M  = 


U1V1 


U?V? 


Ur  Vr 

5  5 


5  J 


One  thus  finds  the  representation 


1  /2 


2 


h(t,q)  =  u  L[l,u,v,uv,v']A 


h5 


where  A 


1 C ( 1 / 2 ) u  1 /2  ,  ( 3/2)u1 72  ,  ( 1 /2)u  1 /2v,(3/2u1 /2V,(1 /2)u  1  v)cos  a 


+  [0,0,u17^,u372,2u172v]sin  a } A 


|-[(1/2)u~1/2,(3/2)u1 /2,(1/2)u  1 /2v, ( 3/2)u1 /2v, ( 1 /2)u  3/2v2]sin  a 


rn  .  1/2  3/2  0  1 /2  n  i ,  . 

+  [0,0, u  ,u  ,2u  v ] cos  a}A  . 


For  the  subsequent  integrations  the  origin  of  the  coordinate 
system  is  shifted  to  the  control  point.  The  new  coordinates  are 
denoted  by  u  and  v. 


U  =  U 1  +  u 


u  =  u  -  u, 


V  =  v1  +  V 


V  =  V  -  V. 


For  each  of  the  six  triangles  a  separate  coordinate  system  (u,v) 
is  introduced  so  that  the  side  !L ,  £  +  1  is  a  line  u  =  const.  For 
1=7  the  point  (H+1)  is  identified  with  point  2. 


We  define  for  a  triangle  1,  l,  1+1;  l  =  2. ..7 


cos  a  = 


v  -  v 
2.  +  1  Z 


r -  ,2  -.2,1/2 


■V  *  <Vi 


'F1' 


*  .  V* 

a-*- 


u  =  u  cos  a  +  v  sin  a 

v=-usina+v  cos  a 


u  =  u  cos  a  -  v  sin  a 

v  =  u  sin  a  +  v  cos  a 


(34) 


w  =  v/u 


(35) 


In  the  triangles  1,2,3;  1,3*4;  and  1,4,5  (Figure  4)  the 
integrand  is  smooth  and  the  same  procedure  as  for  triangles  at 
inner  control  points  is  applicable.  the  region  of  integration  in 
the  uw-system  is  already  a  rectangle.  To  obtain  -1  and  +1  as 
limits  in  both  directions  one  sets 


u  =  (u^/2)  ♦  u(ua/2) 
w  =  ( w  j,  + 1  +  wji)2  +  w(wn  +  i  "  w*)/2 


Then  one  has 


3U,n) 

/V  A 

9 (u , w  ) 


u(1  /4)ua(wjl  +  1 


) 


For  the  triangle  1,6,7  the  transformation  ( Eq .  (34))  gives 

u  =  -  u 
v  =  -v 


As  bef ore , 


w  =  v/  u 


The  limits  for  u  are  0  and  u, ,  and  the  limits  for  w  are  w,  and  w„ 

(M  (?)  6  7 

In  the  representation  for  h  and  h  J  the  factor 
-1/2  -  -  1 /2 

u  =  (u1  -  u)  occurs.  It  is  counteracted  by  the 

transf  ormati on 
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The  limits  for  p  are  u. 


(lower  limit)  and  0  (upper  limit).  To 


transform  in  both  directions  to  limits  - 1  and  +1  we  set 


„  '  1  /2  ,  "  1 /2 
P  =  (u1  ) / 2  -  p(u1  /?.) 


w  =  (  w  +  Wg  ) / 2  +  w ( w  -  w6  )/2 


Then 


“TT'  ■  “,1/2("7  -  V 

3  (p  ,w  ) 

For  the  triangle  1,5,6  one  again  carries  out  the  transforma¬ 
tion  to  u  and  w,  Eqs.  (34)  and  (35),  one  has 

3  (  6 ,  n)  =  - 
3  ( u  ,  w  ) 

The  limits  of  integration  are 


wlower  w5 
ulower  ”  ^ 


W  =  W 

upper  w6 
uupper  u5  U6 


Notice  that  sin  a  <  0  because  u^  <  u^  .  At  point  6,  u  =  0,  but 

u  =  u1  +  u(cos  a  +  w  sin  a) 


The  expressions  for  h^2^  and  h  ^  contain  a  factor  u~1/2,  which 
introduces  a  singularity  at  point  6.  To  obtain  a  smooth  integrand 
we  replace  the  variable  of  integration  u  by  a  new  variable  p 


•' w v. ^ ■>  ■  ■ .»  ".«  *-»  V  ’■’V  *.»  "Jr  P.y  ■>  »  V*3r»  V ■» 


Then 


2 

p  «  u 
2pdp 


+  u ( cos  a  +  w  sin  a) 
=  du(cos  a  +  w  sin  u) 


3U, n) 
3Tp  ,  w  )' 


_ 2pu _ 

cos  a  +  w  sin  a 


The  transformation  fails  if  the  denominator  in  the  Jacobian 
vanishes,  that  is,  for 


w  =  -cot  a 


Then,  according  to  Eq .  (35), 

v  sin  a  +  u  cos  a  =  0 


and  from  Eqs .  (34) 


u  =  0 

This  value  of  u  lies  within  the  region  of  integration.  Let 

w  =  (w6  -  cot  a)/2 

We  divide  the  region  of  integration  over  w  into  two  regions: 

w._  <  w  <  w  and  w.„  <  w  <  w,  .  The  integration  is  carried  out  in 
b  ID  10  b 

the  first  region  in  terms  of  u  and  w;  in  the  second  one  in  terms 
of  p  ana  w. 

In  the  region  w.  n  <  w  <  w,  ,  the  limit  of  p  corresponding  to 
-  1  u  -  D  1/2 

the  lower  limit  of  u,  namely  u  =  0,  is  p  =  u^  .  The  limit  of  p 
corresponding  to  the  upper  limit  u  -  u^  is 

pu pper  =  K  +  u6(cos  a  +  w  sin  a)]  (36) 


It  depends  upon  w. 


0 


=  0  =  u i  +  Ug(cos  a  +  Wg  sin  a) 


Therefore  u^  can  be  replaced  by 


u1  =  u^(cos  a  +  Wg  sin  a) 


Substituting  this  form  of  into  the  expression  (36)  for  the 
upper  limit,  one  obtains 


upper 


-  1/2,- 

u6  ("6 


-  ,1  /2  , 
w)  ( 


sin  a ) 


1  /2 


(We  had  observed  above  that  sin  a  <  0.) 

We  have  generated  a  smooth  integrand,  but  the  boundary  of 
the  region  has  a  singularity  for  w  =  w^.  To  counteract  the 
singularity  in  the  limit  of  w  we  introduce 


c,2  -  i6  - 

w  ,  w  =  Wg  - 

CM 

c r 

2qdq  =  -dw 

Then 

3  ( 5  »  n) 

ilpqu 

3  ( p  ,  q  ) 

cos  a  +  w  sin 

a 

The  limits  for  p 

The  limits  for  q 

obtain  limits  ±1 

1  /2 

are  p,  =  u.  ,  p 

lower  1  _  upp^ij 

are  q,  =  (  wc  -  w,  ~  ) 1  , 

M1  ower  6  10 

in  both  directions  we  set 

=  Ug  q( -  si n 

,  q  =  0 . 

Mupper 

S)’/2 

To 

,  1/oWl  1/2  1/2  ,  .  1  /2 ,  1/2  ,  .  ',1/2 
p  =  ( 1 / 2 ) ( u 1  q(-sin  a)  )  +  p(u^  q(-sin  a)  -  u1 

q  =  (1/2)(Wg  -  w1Q)1/2  -  (1/2)q(wg  -  w1Q)1/2 


T  hen 


jU ,  n? 

Ci  (  r\ 


pqu  r  -  1/2  /  .  ',1/2  1  /?  1,  - 

cos  a  +  w  sin  a  6  M  1  ‘  j 


1  /  2 


1  O' 


The  integration  is  now  carried  out  in  terms  of  p  and  q. 


For  w._  <  w  <  w.^,  one  has  a  rectangle  on  the  u,w-plane.  To 
b  1  0 

obtain  limits  -1  and  +1,  we  set 


u6  ~  ^6 
u  -  2  +  u  2 


W1  0  +  w5  .  ~  W1 0  w5 
w  -  2  +  w  2 


Then 


3U,h) 

3 (u ,w  ) 


¥  u6 (wl 0 


) 


The  integration  is  carried  out  in  terms  of  u  and  w. 

An  analogous  procedure  is  carried  out  in  the  triangle  1,7,2. 

i  /p 

The  critical  term  is  again  u  where, 


u  =  u1  *  u(cos  a  +  w  sin  a) 

the  factor  of  u  vanishes  for  w  =  -  cot  a.  (Here  sin  a  >  0.)  The 
region  of  integration  in  the  w  direction  is  divided.  Let 

wio  °  ^w7  ”  GOt  °^/2 

If  w^  >  w9  (which  is  unlikely),  then  there  is  only  one  region  of 
integration  w^  <  w  <  w^  and  w1  ^  <  w  <  vi^.  In  the  first  region  we 
introduce,  as  before,  instead  of  u 

p^  =  Uj  +  u(cos  a  +  w  sin  a) 

The  limits  for  p  are 

1  /2 

^lower  U1 


pupper  ‘  Cul  *  V003  "  *  “  sin  “]1/2 
-  G7'/2(sin  S),/2<a  -  i/'2 


Next,,  one  removes  the  singularity  in  the  upper  limit: 

2 

q  =  w  -  w„ 


Then  q,  =0  and  q  =  (w...  -  w~) 

Mlower  Mupper  10  7 


A  f urther 


transf  ormation 


(1/2)(w1q  -  w?)1/2  +  q(1/2)(w10  -  w?)1/2 


=  (  1  /2)  [  u, 


,v  r  1/2  -  1  /2,  .  -,i1/2 

’ ) L  u1  +  Uy  (si n  a) j  q 

1/2,  .  1  /2„  1 

p ( 1 /2 ) [ u7  (sin  a)  q-u 


gives  in  both  directions  of  the  p.q-plane  the  limits  are  -1  and  +1 


n)_ 

3  ( p ,  a ) 


pqu  r-  1  /2  ,  .-,1/2  1  /2 1 r 

- - -  [u  (sin  a)  q-u  J  [ ' 


oo s  n  +  w  sin  a 


1  0  T 


In  the  region  Q  <  w  <  w^,  we  introduce  u  and  w  by 

u  =  (1  / 2 ) u  ?  +  u( 1 /2)Uy 
w  =  ( 1 / 2 ) ( w ^  +  w(1/2)(w2  -  w  ) 

In  the  u,w-plane  the  limits  are  -1  and  +1  in  both  directions.  One 


3  ( £  ,  n ) 
3  ( u  ,  q  ) 


=  (u/4)u7(w0  -  w  ) 


For  control  points  immediately  upstream  of  the  trailing  edge 
f i-i^-1 ) ,  one  determines,  as  for  other  hexagons 


Xm(  J, )  ,  n  ( {, ) 


C.  =  1  ...  7 


n.  =  (y_ 

m(n),n(0 

Next,  one  shifts  the  coordinate  system  to  point  3 

n£  =  nt-n3 

The  sweep  angle  at  the  trailing  edge  is  given  by 

-  2  -2,1/2  .  -  ,  ,  -  2  -2,1/2 
cos  a  =  n^/U^  +  )  ,  sin  a  =  ^/U^  +  ) 

A  coordinate  system  is  introduced  so  that  the  trailing  edge  is  the 
line  u  =  0 

\  =  u  cos  a  +  v  sin  a  ,  u  =  £  cos  a  -  n  sin  a 

n  =  -u  sin  a  +  v  cos  a  v  =  £  sin  a  +  n  cos  a 

Next  one  determines  the  values  of  u^  ,  v  ,  2.  =  1...7, 

(u^  =  0 ,  u  =  0).  The  steady  problem  suggests  a  representation  for 
h  of  the  form 

h(£*n,n)  =  [i,u,v,  uv,  v2,(-u)^/2] 

For  the  reason  given  in  conjunction  with  the  presentation  of  h  in 

2 

the  vicinity  of  leading  edge,  the  term  u  is  not  included  in  the 
row  vector  on  the  left. 

In  the  present  form,  only  6  functions  of  u  and  v  and  6 
coefficients  a^  are  available;  therefore,  only  an  approximate 
matching  is  possible.  The  procedure  is  analogous  to  that  for  an 
inner  control  point.  One  introduces  a  7  by  6  matrix 

69 


»  ,  ■  _  «  ,  • 


and  f orms 


A  =  [M  +  M]  ^ 


One  thus  obtains  the  approximation 


h ( t  ,  C  .  n) 


[1 i u, v ,uv. 


v 


,  3/2  , 
u) J  ] 


■h1  (  t  )" 

hy  (  T  ). 


Hence , 


h ( 2 } ( t , £ , n)  =  {[0, 1 ,0,v,01”3/2(-u)1 /2]cos  a 


+  [ 0 , 0 , 1 , u  ,  2v  ,  0  ]  s  i  n  a]}  A 


■h1  (x  n 


L  h  ^  (  T  )_ 


and 


h ( 3  ( r.C.n)  =  {-[O.I.O.v.O  -(3/2) ( - u  )  1 72 ]  sin  a 


+ 


‘  h  1  (  t  )“ 
_  h ..  (  t  )_ 


[  0 , 0 , 1 , u , ?v , 0]cosa] } 


A 


For  the  integration  the  origin  of  the  coordinate  is  transferred  to 
the  point  1,  (notice  that  u.,<0) 


u  =  u  -  u i  u  =  u  ^  +  u 


v  =  v  ~  v  1  v  =  v1+v 


The  six  triangles  1,2.,  2  +  1  are  treated  separately.  In  each 


triangle  a  coordinate 

system 

is  introduced 

so  that 

U2 

U2+1 

cos  a  =  (ui+1- 

V 

/ 

1 

+ 

1  =5 
' — «' 

i — i 

V2  * 

1 

+ 

1  > 

’V2 

)2:1/2 

s  i  n  r.  =  (  v  t  + 1  - 

V 

/ 

[<ui.r 

u(F  • 

( Vr 

'V2 

)2],/2 

u  -  ucosa+vsina 


u  =  u  cos  a  -  v  sin  a 


v  -u  sin  a  +  v  cos  a 


v  «  u  sin  a  +  v  cos  a 


One  determines  ug  =  u&+1,  v£ ,  and  vJ+1. 
Next  one  introduces 


w  =  v/  u 


and  determines  w.  and  w.  , 

2  2  +  1 


Then 


3  (  6 »  n)  =  - 

9  ( u  ,  w  ) 


In  the  triangles  156, 
region  of  integration 
rectangle.  To  obtain 
coordinates  u  and  w; 


1 6 7 ,  and  172  the  integrand  is  regular.  The 
for  these  triangles  in  the  uw  plane  is  a 
limits  -1  and  +1,  one  introduces  new 
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+  w 


+  w 


w  = 


( w  +  w„  ,  ) 
2.  1  +  1 


2+1 


i  nen 


3U,_n) 
3 (u ,w  ) 


u 

11 


“t(wt.r  "2 


-  -  ( 2 )  n ) 

For  the  triangle  13^.  u  =  u,  v  =  v.  In  h  and  h  ,  a  term 
1  /? 

(-u)  “  occurs,  which  introduces  a  singularity  into  the  integrand 

One  has 


Therefore,  we  set. 


-  u  =  -  ( u1  +  u )  . 

p2  =  -  ( u1  +  u ) 
2pdp  -•  -  du 


Then 


3  (  £  ,  n) 

3  Cp  ,  w ) 


-2up . 


The  limits  for  p  are 

,  >1/2 
p ,  =  (  -  u  .  ) 

p lower  1 

p  =  0 . 

Ku  ppe  r 

Limits  ± 1  are  obtained  by  introducing 

p  =  ( 1  /2) (-u1  )1  /2-  ( p/2) ( -u1 ) 1  12 

w  =  1/2  (  w  +  w  ,  )  +  (w/2)(w,-w,). 


T  hen 


3  ( £ , n)  =  up 
3(p,w)  2 


(  -u. 


w. 


) 


For  the  triangle  1,  2,  3.  one  has  after  the  initial  transformation 
(-u)  =  -U.J-U  =  -u^-  u  (cos  a  +  w  sin  a) 
here  sin  a  >  0,  since  >  u9. 

1  /2 

To  counteract  the  singularity  introduced  by  (~u)  ,  we  set 

2  - 

p'  =  (-u)  -  - u ^  -  u  (cos  a  +  w  sin  a) 

Here  sin  a  >  0  because  u^>u^-  The  transformation  fails  for 

w  =  -  cot  a. 


Let 


w1Q  =  (  w  ^  -  cot  a )  / 2  . 

If  w^>  w  ,  there  is  only  one  region  of  integration 

If  w9  <  w  ,  one  has  the 

<  w 

For  w <  w  <  w  q  ,  the  integrand  is  smooth  and  the  integration  is 
carried  out  in  terms  of  u  and  w.  To  obtain  limits  O,  we 
introduce  instead  of  u  and  w. 


w0  <  w  <  w 

^  5 


two  regions 


<  w  and  w  <  w  <  w 
I  J  1  U  j 


-J  -J  ^  -J  in '  U"-1  wir  r.i’.-y.'  wj  TV  wj  ■  u.»  v  y.’vv.Tjr.v  '.■*  7« v»  t  '  ” ' "  ■  ■•  tt y T^rysr -,v .t- 


u  =  u(u2/2)  +  u(u,>/2) 


w  =  (w1Q+  w2)/2  +  w(w1Q-  w2)/2) 


Then ,  for  w,  <  w  <  w 


1  0 


3Uj.rO 

3 (u  ,w  ) 


“o(w,  «-wJ 


IT  u2  ^  1  0  2 


For  w1 q  <  w  <  ,  p  is  introduced  by  the  above  formula.  As 

limits  of  integration  for  p,  one  finds 


^lower  ~  ^  u 1 ^ 
P 


1  /  2 


upper 


( -  u1  -  u^(cos  a  +  w  sin  a)) 


1  /  2 


Rut  since  u  =  0  for  u  =  ,  w  =  ,  one  has  -u1 =  u^,  and  then  from 


the  transformation  formulae 


-u1 =  u^(cos  a  +  w3  sin  a) 


Theref  ore , 


„  /  ~1  /2  ,  .  . 1  /2  -.1/2 
pupper  ‘  (u3  (31n  a))  (W3  '  w) 


Let 


-  w  =  q 


‘lower 


( w3  wi o  ^ 


1  12 


for  wi0  >  w2 


(W  -  w2) 


1  /2 


for  w1  <  w 


upper  »  0 
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V, 

•1 


*3 


?3 


vJ 


To  obtain  limits  +1,  we  introduce  new  variables  p  and  q  by 


,  .  -.1/2  ,  .1/2.  ,  .  1/2,  .  -.1/2 

(sin  a)  q  +  C-^)  ]  +  (  1  /2)p[u^  (sin  a)  q 


.1/2  ,  .  ~  ,  -  -  .1/2 

w-  -  w13)  -  (1/2)q  (w^  -  w1 ^ ) 


nen  , 


3_Ui_ni 

3(p,q) 


__  ._u._P._q__ 


cos  a  +  w  sin  a 


r  ,  -.1/2  -1/2,  .  -.1/2  nr  - 
=  [(-u)  - u  ^  (sin  a)  q  ] r  w  ^ ■ 


If  w?  >  w^,  one  replaces  in  these  formulae  w^  q  by  u 


For  the  triangle  1,  4,  5,  one  introduces 


w  «  (w^  -  cot  a)/2 


and  obtains  regions  of  integration 


w4  <  w  <  w5 


for  w1Q  >  w^ 


and 


w,  <  w  <  w  arid  w.  <  w  <  wc  t'or  w,  <.  w, 
4  ID  10  p  1  J  , 

For  *  he  region  w  <  w  <  wr  (which  occurs  in  the  so  rond 
o  d t a i n s  )  i m  its  +1  by  introducing  n ew  coo r d i n a t es  u ,  w 


p--» ■  »  -  «  -  vi  ■  ■■ VTV v»Y-« v»  '.~v  V« V~v yv ».^1V'\.y'.V'^^.,S 


w2  =  (w5+  w10)/2  +  w(w5-  w10)/2 


3 

.1 

."i 


One  has 


3  ( S  ,  n)  _  u 


3(u,w) 


ur(wc  -w1n) 


IT  u5 v  5  "W10 


.* 

•/ 

V 

V 

u 

M 


In  the  regions 


One  sets 


w4  <  w  <  w10 


(  -  u .. ) - u  (eos  a  +  w  sin  a) 


Here  sin  a  <  0.  The  limits  of  p  are 


i  ower 


( -u  1  ) 


1  /  2 


r,  -  ,  -  -  ,  1  /2  -1/2,  .  - .  1  /2  ,  -  -  ,1/2 
pupper  =  u4(cos  a  +  w  sin  a)  J  =  (-sin  a)  (w-w^) 


The  form  of  p  suggests  that  one  introduce 

upper 


q  =  (w  -  wn) 


t  h  e  n 


q.  =0 
dower 
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# 
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.  ^ 

>S 
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v 

« 


f or  w 


q  is  (w._—  w,.) 

Mupper  10  4 


1  0  <  W5 


,  -  -  ,  1  /2  _  -  .  - 
(w1Q-  w5)  for  w10  >  w5 


One  introduces  new  variables  p  and  q  to  obtain  limits  ±1 
p  =  (1/2)[uli1/2(-sin  a)1/2q  +  (-  u1  )1//2]  +  (1/2)p[u^1/2(-sin  a)1/2q-(-u1  )  1  /2] 

q  =  (1/2)([w10  -  w/72  +  (1/2)q  [w1Q-  w^)1/2] 


then 

3 ( 5, n) 

- - « —  — 

9 ( p  ,  q  )  cos  a  +  w  sin  a 


u  p  q 


r,  -  J/2  -1/2,  .  -  N 1  /2  ir~  -  ,  ,1  /2 

[(-u^  -u^  (-sin  a)  qltw^w^)] 


If  w1Q  <  w^,  one  applies  the  above  formulae  with  w1Q  replaced  by 


w5. 


For  a  control  point  at  the  trailing  edge  (i  =  i t )  , 

(Figure  12),  one  first  determines 

^  l  ~  X  m  ( l )  ,  n  (  H ) 

\  =  ymU  )  ,  n  (  i  ) 

In  introducing  a  triangulation  of  the  wake,  it  facilitates  the 
determination  of  h  if  one  lets  n-^r^.  n4  =  n1  .  This  means  that 

the  basic  quadrangles  are  trapezoids.  Next  we  shift  the 
.'card  1  cat-/  systems  into  point  1 

7  7 


One  introduces  a  Cartesian  system  u,  v,  where  u  is  perpendicular 
to  the  trailing  edge 

-  I/O 

cos  a  =  (n5-n2)/U5-S2)  +(n5-n2)  ] 

sin  a  =  (^5-C2)/[C5-C2)2+(n5-n2)2]1 /2 

\  =  u  cos  a  +  v  sin  a  ,  u  =  £  cos  a  -  n  sin  a 

n  =  -u  sin  a  +  v  cos  a  ,  v  =  \  sin  a  +  n  cos  a 

The  steady  problem  suggests  a  representation 

2 

h  =  [1,  u,  v,  uv,  v  ,  0] 

2  1/2 

for  the  wake,  and  h  =  [1,  u,  v,  uv ,  v  (-u)  ]  for  the  wing. 

2 

A  term  u  has  not  been  included  for  reasons  given  above. 
Again,  only  approximate  matching  is  possible.  The  procedure  is 
analogous  to  that  for  an  inner  control  point.  One  substitutes  u 
and  v  for  points  1  to  7  in  the  expressions  applicable  either  for 
the  wings  or  the  wake,  and  obtains  a  7  by  six  matrix  M 


pnnm  *"y  Ti  v"j 
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M  = 


1 

0 

0 

0 

0 

0 

1 

0 

V2 

0 

2 

V2 

0 

1 

u3 

V3 

U3  V3 

2 

v3 

0 

1 

u4 

v4 

U4V4 

2 

v4 

0 

1 

U5 

V5 

U5V5 

2 

V5 

0 

1 

u6 

V6 

U6V6 

2 

v6 

,  ,  3/2 

”  u6 

1 

U7 

v7 

U7V7 

2 

V7 

,  ,3/2 

(  -  u,. )  J 

Then  one  f orms 


-1 


A  =  [M+M]  M+ 

and  obtains  the  approximation 


a  1  (  t  )~ 

=  A 

h^  ( -t  ) 

_a7(T)_ 

hrj  (  T  ) 

Then ,  in  the  wake 


i/2^=  {  [  0  ,  1  ..  0  ,  v  ,0 , 0]cos  a  +  [ 0 , 0 , 1 , u  ,  2v  ,  0]si  n  a]}A 


h1  (  t  ) 


h  y  (  t  ) 


and,  at  the  wing 
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il 


,  r 
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; 


1  /2 


hi,  =1  CO,  1  ,  0,  v  ,0,  -  (3/2)  (-u)  1  ' 11  ]cos  a+[ 0 , 0 , 1  ,u,2v,0]sin  a]}A 


h^  ( T 


Furthermore,  in  the  wake 


hp^  =  { -[0 ,  1  ,  0 ,  v  ,  0 , 0]si  n  a  +[  0 , 0  ,  1  ,  u  ,  2v  ,  0  ]  cos  a]}  A 


h  i  ( t ) 


hy  (  T  ) 


and  at  the  wing 


V 

h^  =  {  -  [0  ,  1  ,  0  ,  v  ,  0  ,  -  (  3/2  )  (  -  u) 1 /2  ]s  i  n  a+[ 0 , 0 , 1 , u , 2v , 0 ] cos  a]}A  : 

h  rj  ( 

Next  one  considers  individual  triangles  and  rotates  the  uv  system 
into  uv  systems  so  that  the  sides  i,  {,  +  1  are  lines  u  =  const. 


cos  a  = 

(\.i  -  Vlt'Vr 

V2* 

(v».r 

V2 

],/2 

sin  a  = 

(u,.i  "  VK'V, ' 

V2* 

<vt.r 

v2 

],/2 

=  u  cos 

a  +  v  si n  a  , 

u  =  u 

cos  a 

-  v  s  i  n  a 

v  =  -usina+v  cos  a  ,  v=usina+v  cos  a 


~  V 

Moreover,  w  =  —  . 

w 

In  the  triangles  123,  13^,  and  1^5,  the  integrands  are 
regular.  To  obtain  limits  of  integration  ±1,  one  introduces  new 

A  A 

variables  u  and  v  by 


u  =  ( u^/2 )  ■*  u ( u^/2 ) ) 

w  =  (w^+1  +  w^/2  ♦  w  (w£  +  1  -  we)/2 

One  has 


7"^  ■  ozo  G  -  “t> 

3 (u  ,w  ) 

In  the  triangles  (156),  (1,6,7),  and  (1,7,2), 

1  /? 

representations  of  h  includes  a  term  (-u)  ',  which  introduces  a 

singularity  long  the  line  2,  1,  5.  One  has 

(-u)  =  -u(cos  a  +  w  sin  a). 

In  the  triangle  (1,6,7),  u  is  negative.  Therefore, 
cos  a+wsina<0 

,  ,1/2  -  1  /2r  ,  ~  -  .  n1/2 
(-u)  =  u  [-(cos  a  +  w  sin  a)  1 

2  ~ 

We  set  p  =  u  ,  and  leave  w  unchanged.  To  obtain  limits  ±1,  a 
p,  w  system  is  introduced  by 

-1/2  ~  ~  i  /? 

p  =  (u6)  ~/2  ♦  p(u6) 1 'd/2 

w  =  (w r  +  wr.)/2  +  w  (w,-  w._)/2) 

6  5  6  5 

t  hen 

3U,h) 

- K - - 

3 (p , w  ) 

In  the  triangle  (1,5,6) 

u  =  0 ,  for  w  =  w  ^ . 


up  ,  -  ,1/2,- 
2  u6  w6 


-  w5) 
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-V.‘ 


V  s 


Therefore,  cos  a  +  w,sin  a  -  0 ,  and  ( -  u )  =  -  u  (  ( w  -  wr)sin  a.  Then 

j  J 

one  introduces 

o  ~ 

Q  “  =  (  W  -  W  r  ) 

and  subsequently,  q  and  u  by 

q  =  ((w^-  )  /  2)  *  q(w^  "  wt-)  /2 

u  =  (Up/2)  +  u  (u,  )/2. 
o  b 

Then , 

=  uq  (1/2)Uc  )[vA-  vi^)W2 

3(p.q)  5  6  ’ 

In  the  triangle  1,7,2,  one  sets 

2 

q  =  ( w  0  -  w ) 

and  introduces 

q 

u 

Then 

~ —  =  uq(1/?)u  (w  -  w)1^ 

3(u,q)  1  '  1 

For  a  control  point  at  the  x  axis  (the  axis  of  symmetry), 
one  has  in  principle,  the  same  procedure  as  for  an  inner  control 
point,  but  it  is  possible  to  take  the  symmetry  properties  of  the 
flow  field  into  account.  Referring  to  Figure  1  3,  one  has 


I/O  -  -  .1/0 

=  (  (  w  0  -  w  )  /  2)  ~  q  (Wj  -  w  ^  “ / 2 

=  u  /2  +  u ( uy/2 ) 


m.|  =  i  =  0 

m^=i  +  1  n  ^  =  0 


mc  =  l 
5 


n5  =  1 


m6  =  i_1 


n6  =  1 


-  i-1 


=  0 


=  x 


m ( SL )  ,5(1) 


-  ^  -  y 


mil) ,n(L) 


*•  -  1  ,  4,  5,  6, 


Moreover 


The  origin  of  a  £,n  system  is  placed  into  point  1 


For  a  flow  field  which  is  symmetric  with  respect  to  the  x 
one  sets 


h  =  [1  ,  l,  l2 ,  n2] 


To  match  at  five  points  1,  4,  5,  6,  and  7,  one  introduces 


axis, 


a  matrix 
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and  forms 


A  =  [ M  +  M]  1 M + 


Then  one  obtains  the  following  approximation  for  h 


h-  [i ,l,l2 ,n2) 


r!jl 

k 

h7 


Hence, 


h(2) 

symmetric 


[0, 


2£,0]  A 


N  1 
N 

h5 

*"  i 


h 


(3) 

symmetric 


[0,0,0, 2n]  A 


L1J 


For  the  antisymmetric  part  of  h,  one  sets 


h  .  .  =  [  u  ,  £  n] 
anti 


and  definei 


n^  »  ^5^5 


n6  •  56n5 


A  =  M 


now  , 


anti 


[n,  5n]  A 


anti 


[ 0 , n]  A 


anti 


[0,53  A 


The  steps  for  the  individual  triangles  are  the  same  as  for  the 
interior  of  the  wing.  Only  the  triangles  in  one  half  of  the  wing 
need  to  be  considered.  The  sign  of  the  contribution  of  the  other 
triangle  are  determined  by  the  symmetry  cr  antisymmetry. 


As  mentioned  in  Section  IV,  the  leading  edge  has  been 
modified  in  the  vicinity  of  the  axis  of  symmetry  (see  Figure  2) 
For  a  control  point  next  to  the  leading  edge,  one  obtains  a 
configuration  consisting  of  5  triangles  (Figure  I1!). 


n1  =  1 


n  1  =0 


m4  *  2 


n4  =  0 


-  ■'  s  , <■  i *-«y.  r.  <■.  j-m  s,  .'.v'AYi/.v.v.v'w;.- . 


mc  *  1 

j 


n5  -  1 


m6  =  0 


n6  =  1 


^  -  x(rv  V* 


1^=  y  (m  ?,  * n  2, )  -  1  =  1  .4,5,6 


^2  =  ^6  ’ 


V  "n6 


^3  *  ^5 


V  -n5 


We  introduce  an  auxiliary  point  6' 

^6 '  =  ^6  ’  n6 '  =  0 
and  place  the  origin  into  point  6' 

^ l  ~  ^  l  ~  ^6 '  '  r|2.  =  riS, 

For  the  symmetric  case,  h  is  then  represented  by 


h  ,-r  1/2  r-3/2  ;i/2  2,  a  ' 
n  =  ( C  ,  S  n)  2 

a3 


then  hg  =  0. 


3  1,4, 

and  5, 

i;/2 

v3/2 

v/2 

r  3  /2 
^4 

iy 2 

H'2 

=  1/2  2 

55  n{. 


and  evaluates 


A  =  M 
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-  ^  n."  v*  v.'.  -Ai\' 


k-  _v 

.'.V 

V.y'Vl 


symme t r i c 


'  L  s  1  •?  »  S  1  ’  J  4  I  “  JJ 

K 

For  the  antisymmetric  case,  one  has 


anti 


r  1  /  2  ,-1/2 

=  5  n(c5  n5. 


hence  ,  h 


12) 


and  h 


(3) 


After  the  introduction  of  the  point  6',  one  deals  again  wit 
a  hexagon  in  which  the  point  6'  plays  the  role  of  point  7.  The 
only  difference  lies  in  the  fact  that  the  line  6'2,  as  well  as 
6,6’,  are  part  of  the  leading  edge.  The  treatment  of  the 
triangles  (1,4,5),  (1,4,6),  and  (1,6,6')  is  the  same  as  for  a 
control  point  adjacent  to  the  leading  edge.  The  contributions  of 
tne  remaining  triangles  are  found  immediately  on  the  basis  of 
symmetry  consideration. 

The  treatment  of  the  trailing  edge  is  different,  because 
there  exists  control  points  directly  at  the  trailing  edge.  On  th 
other  hand,  the  singularity  at  the  wing  has  the  character  of  u^' 
while  the  flow  in  the  wake  is  free  of  a  singularity. 

The  author  proposes  not  to  modify  the  trailing  edge  at  the 

3/2 

line  of  symmetry,  and  to  disregard  the  u  "  singularity.  For  the 
control  points  at  the  line  of  symmetry,  next  to  the  trailing  edge 
as  well  as  at  the  trailing  edge,  h  is  then  represented  by 


[1 


r  r2  r  21 

,  6  ,  n ,  6  ,  6  n ,  n  J 


i.e.,  in  the  same  manner  as  any  other  control  point  in  the 
interior  or  the  wake.  The  symmetry  relations  can  be  taken  into 
account  as  for  other  control  points  on  the  x-axis.  One  might 
experiment  with  other  approximations,  but  it  is  likely  that  the 


reatment  of  a  root  region  will  always  remain  somewhat 
■a  ns at i sf ct or y  until  one  extends  the  method  to  include  the 
and  the  fairings  at  the  root  sections.  This  was  already 
mentioned  in  Section  IV. 


SECTION  X 


THE  NUMBER  OF  PIVOTAL  POINTS  TO  BE  USED  FOR  THE 
INTEGRATION  WITHIN  ONE  TRIANGLE 

The  number  of  pivotal  points  to  be  used  for  the  integration 
depends,  of  course,  on  the  variations  of  the  integrand  within 
this  region.  These  variations  are  brought  about  mainly  by  the 
denominator 


DEN  =  p[MU-x)  ♦  p)2 


Here 


p  =  [ ( £-x .  .  )2  +  ( n  -  yi j )  ]' /<~ 

and  x  =  x i j ,  y  =  y  —  are  the  coordinates  of  the  control  point 

under  consideration,  £  and  n  are  the  umbral  variables  of 
integration . 

It  is  proposed  to  take  the  ratio  of  the  maximum  to  the 
minimum  over  the  triangle  of  the  expression  DEN  as  criterion  for 
the  choice  of  the  integration  formula.  The  extrema  occur  at  the 
contour,  but  not  necessarily  at  the  corners  of  the  triangle.  We, 
therefore,  compute  DEN  at  the  corners  of  the  triangle  ard  at  the 
midpoints  of  the  sides,  approximate  DEN  along  each  side  by  a 
second  degree  polynomial,  and  then  determine,  separately  for  each 
side,  the  extrema.  The  corners  are  assigned  indices  1,  2,  and  3 
and  the  midpoints  of  the  sides  from  l  to  i  +1  the  index  t  +  1/2. 
'For  i  =  3 1  we  identify  8.  ♦  1  with  1  .  ) 

The  following  steps  are  then  carried  out.  One  determines 
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■V 
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^^7  ’  yyy  v”* ;.  ^  .^>  •;*?  > ,  a  • 


P: 


^2  +  1  / 2  =  (  1  /2)  (^2  +  ^ 8.  + 1  ) 

ni-1/2  =  C 1 / 2 ) ( n  ^  +  yi  +  1  ) 
for  2.  =  1  ,  2 ,  3 , 


<-  i  ^ 


and  evaluates  D  E  N  ^  and  D  E  N  ^  +  -|  ^  2 

Considering  a  side  (2., 1  +  1),  we  introduce  an  auxiliary 
coordinate  z  so  that  z  =  -1,  z  =  0,  and  z  =  +1,  respectively,  at 
the  points  2, 2  +  0/2)  and  2,  + 1  .  Then  one  has  an  approximation  for 
DEN  along  the  side  (2,2+1) 


DEN(z)  =  DENj,  +1  /2  +  a  1  z  +  a2z‘ 


with 


a1 *  (DEN^  +  1 -  DEN ^ ) / 2 


a2=  (DENt  +  1-  DENr  2DEN p_  +  1 /2  ) /2 


The  extremum  lies  at 


z  = 


a1 

2a. 


If  a  =  0,  it  is  practical  to  replace  it  by  a  on 

10°  say,  to  avoid  ar;  exceptional  treatment  ‘ 
The  extremum  lies  within  the  side  of  the  tri 


|  z  |  <  i  . 


It  is  given  b- 
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ext , i , fc+1 


DEN 


If  |z|  >  1 ,  we  set  DENext  t  l+1  •  DEN^.  Then  one  has  for  the 
overall  maximum  and  minimum 

DEN,nax  ‘  ”ax<DENf  DENext  =  '  '  '  •  2’  3) 

DENmin  ‘  "ln<DENf  DENext  1  '  ’•  2’  3) 

The  ratio  DEN„/DENm .  is  used  as  criterion  for  the  choice 
max  min 

of  the  integration  formula.  Of  course,  for  triangles  at  some 
distance  from  the  control  point,  one  can  dispense  with  this 
detailed  evaluation. 


SECTION  XI 


SPECIAL  CHOICE  OF  THE  TIME  DEPENDENCE  OF  THE  UPWASH 

With  the  procedure  described  in  the  preceding  sections,  one 
can  determine  the  function  h  and  with  it  the  potential  4>  and  the 
pressures  generated  by  an  arbitrarily  given  wing  displacement 
g(x,y,t). 

In  aeroelastic  applications,  the  displacements  of  the  wing 
surface  are  given  by  a  moderately  large  number  of  known  functions 
of  x  and  y  multiplied  by  unknown  coefficients  that  depend  upon  t 

n 

g(t,x,y)  -  £  an(t)gn(x,y)  (37) 

In  this  section  the  procedures  will  be  discussed  for 
displacements  expressed  in  this  form.  During  this  discussion, 
the  underlying  integro-dif f erential  equation  is  summarized  in  the 
form 

rlh  +  -► 

SJ  ♦  /  A(  t  )h  (t-x  )dx  =  w ( t )  (38) 

dt  o 

Setting  t  »  t  -  t,  one  obtains,  alternatively. 

♦  J  A(t-T)ii(T)dT  -  w  ( t )  (39) 

dt 

o 

This  equation  is  closely  related  to  Eq .  (7).  The  elements  of  the 

vectors  h  and  w  are,  respectively,  the  time  dependent  parameters 

by  which  h(t,x,y)  (and  with  it  the  potential)  and  the  time 

dependent  values  of  w(t,x,y)  at  the  control  points  are 

characterized.  The  matrix  A  has  elements  which  depend  (in  Eq . 

38)  upon  the  retardation  t.  It  arises  if  one  discretizes  with 

respect  to  the  space  coordinates  x  and  y,  but  considers  x  as 

continuous.  Th's  was  our  original  concept.  In  expressing  the 

coefficient  C,..x  ,  ,  it  had  been  assumed  that  the  matrix 

(  i j ) , (mn  )k 


elements  are  given  in  tables  with  x  as  independent  variable  and 
the  necessary  interpolations  have  been  carried  out.  In  contrast 
A(t)  is  considered  as  a  continuous  function.  One  remembers  that 
the  points  on  the  wing  surface  are  denoted  by  pairs  of 
subscripts.  (This  facilitates  the  identification  of  their 
neighbors.)  If  one  thinks  of  a  matrix  in  C,  one  has  (i,j)  as  row 
and  (m,n)  as  column  indices.  If  one  writes  the  equation  in  the 
above  form,  it  is  assumed  that  one  first  carries  out  the  spatial 
discretization  and  afterwards  performs  the  integrations  over  £ 
and  n  required  in  B^,  B0,  or  Bg  at  the  values  of  x  evaluated  for 
the  pertinent  values  of  (x-£)  and  (y-n).  In  the  procedures 
proposed  in  the  first  part  of  this  report,  the  integration  had 
been  carried  out  over  areas;  the  fact  that  within  such  an  area,  t 
varies  had  been  taken  into  account  at  the  same  time. 

For  displacements  in  the  form  Eq  (37),  the  function 
w(t,x,y),  which  enters  the  aerodynamic  part  of  the  problem  is 
given  by 

n  .  -1  9gn(x,y) 

w(t,x,y)  -  Z(an(t)M  gn(x,y)  +an(t)  - ^ -  )  (HO) 

Here  a’  denotes  the  derivative  of  a  with  respect  to  its  argument. 
In  the  two  dimensional  case,  one  has,  for  instance,  for  a 
plunging  motion 

g  ( t ,  x  )  -  a  ( t ) , 

then 

w(t,x)  -  M_1a' (t) 

For  a  rotating  motion,  one  has 

g(t,x)  -  a(t)x 

then 
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In  the  further  discussions  we  write 


w(t,x,y)  =  1  c  (t)w  (x,y)  (Ml) 

i  +  1  1  * 

In  general,  one  function  g(x,y)  will  generate  two  functions  w£, 
namely  g(x,y),  and  g  (x,y).  The  vector  w(t)  is  formed  from  the 
values  of  w(t,x,y)  at  the  control  points 

w(t )  -  Z  c  ( t )w 
l  1  1 

The  terms  of  the  sum  are  now  considered  separately. 

For  this  purpose,  the  function  c^  is  written  in  the  form 

t+a 

c  ( t )  -  f  c: (o)H(t-o)do  ♦  C . (0)H( t )  ,  a  >  0  (M2) 

1  0  %  * 

where  the  Heavyside  step  function  H(t-t)  is  defined  by 

H  (t-T)  -  0  t  -  t  <  0  (M3) 

H  (t-x)  -  1  t  -  t  >  0 

The  second  term  on  the  right  of  Eq .  (M2)  can  be  omitted  if  one 
postulates  <^(0)  -  0  and  admits  a  discontinuity  of  c^  at 
t  -  0;  c^  will  then  have  a  delta  function  at  o  -  0.  An 
alternative  form  of  c^  (not  used  here)  is 

t+a 

c^t)  »  j  c^  (t )  6  ( t-x  )dT  ,  a  >  0 

It  is  connected  to  Eq .  (M2)  by  an  integration  by  parts.  Notice 
that 


Eq.  (42)  suggests  that  we  wrtt^ 

t  +  a 

i^U)  -  J  c*(o)hi(t,o)do  (45) 

o 

Then  one  obtains  from  Eq .  (39),  by  inspection, 

dh  ( t ,  o )  r 

— -  +  J  A( t-T  )hj (t  ,  o )dx  -  wfc  H(t-o)  (46) 

o 

The  discussions  which  end  with  Eq.  (52)  show  how  Eq.  (46)  is 
obtained  formally.  One  substitutes  Eq.  (45)  and  Eq .  (42)  into 
Eq.  (39)  and  interchanges  the  sequence  of  integrations  in  the 
second  term.  Then  the  integration  over  o  is  the  outer 
integration  in  all  terms,  always  with  a  factor  c^(o).  Since 
c^(o)  is  arbitrary,  the  other  factors  combined  must  vanish. 

We  set 

t  -  o  -  t  (47) 

T  -  0  -  T 

and 

a 

h^U.o)  -  ht(t>o,o)  -  ht(t,o)  (48) 


Then  one  obtains 


t  <  0 


(49b) 


dh. 


dt 


( t ,  o  )  ♦ 


t 

/  A( t-x  )h^ (t  ,  o)di 
o 


0, 


From  the  initial  condition 


ht(0,o)  -  0 


which  becomes 


h^ (-o , o)  «  0 


and  Eq.  (49b),  one  obtains 


hl(0,o)  -  0  (50) 

The  parameter  o  does  not  occur  explicitly  in  Eq .  (49a)  and  in 

A 

initial  condition  (50).  Therefore,  h  does  not  depend  upon  o 

A  A 

hjt.o)  -  hjt)  (51  ) 

Eqs.  (49a)  and  (50)  are  therefore  rewritten 


dh  r 

-=*■  (t)  ♦  J  A(t-T)hl(t)dT  -  ,  1-^(0)  -  0  (52) 

dt  o 

This  equation  is  solved  within  the  desired  region  of  values  t  for 
all  vectors  w^ ,  which  occur  in  Eq.  (41).  Then 

A  A 

hjl(t  ,o)  -  h£(t-o)  (53) 

The  operator  on  the  left  of  Eq .  (39)  is  denoted  by  L 


L(fi)  -  ft  +  1  A(t-T)tl(T)dt 


Eq.  (46)  written  in  the  form 


L(h£(t)  -  ca,  ( t ) 


is  solved  by 


h£(t)  -  L"1 (ct(t)wt 


o’  (o)h. (t-o)do 


(provided  that  0^(0)  -0).  As  was  mentioned  above,  because  of 
this  condition  one  must  allow  that  in  general  c^(o)  has  a  delta 


function  at  o-O. 


Applying  an  integration  by  part  to  the  right  hand  side  of 
the  last  equation,  one  obtains 


ht(t)  -  L  1  (c£(t)wjl) 


|  c. (o)h' (t-o)do 
0  1  *• 


(because  c^(o)  vanishes  from  o-O  and  hfc  vanishes  for  o  -  t) 


The  function  h^  arises  as  one  integrates  Eq.  (52).  Then 


h(t)  -  L~ 1 ( w ( t )  -  I  h.(t)  -  l  j  c. ( o )h I (t-o)do  (56) 

l  1  i  o  x  1 


The  potential  at  the  control  points,  again  written  to  form  a 
vector,  is  then  given  by 


$(  t )  -  2itii(t ) 


The  pressures,  in  dimensionless  form,  are  given  by 


p(x,y,t)  =  ♦(M~1<f>t+  <frx) 


(58a) 


.'tvTvivTv”**1 


where  the  upper  and  lower  signs  refer,  respectively,  to  the  upper 
and  lower  wing  surface.  It  follows  from  Eqs.  (52)  and  (50)  that 

h£(0)  -  w£  (58b) 

The  pressures  enter  into  the  equations  of  motions  for  the 
wing  in  the  form 

I  *  p(t,x,y)qn(x,y)dxdy  (59) 

The  functions  qn  are  identical  with,  or  at  least  closely  related 
to  the  function  gn(x,y)  of  Eq.  (37).  They  are  usually  available 
in  the  analytic  form  so  that  3gn/3x  is  readily  formed.  Using  Eq . 
(58a),  one  then  obtains  for  these  integrals  I 

1  “  M_1  If  4>tqn(x'y)dxdy  +  II  *xqn  (x  , y  )dxdy 

In  the  second  term  an  integration  by  parts  with  respect  to  x  is 
carried  out  (because  (3qn/3x  is  more  readily  available  than  <t>x). 
Taking  into  account  that  4>  =  0  at  the  leading  edge,  one  obtains 

1  '  M~'  H  -  /^trailing  Vxt rail  lng(y  >  )d* 

-  ||  <)»(x,y)qn^x(x,y)dxdy  (60) 

Assume  that  functions  u(x,y)  and  q(x,y)  are  characterized  by 
their  values  at  pivotal  points,  and  thus  result  in  vectors  u  and 
q.  Then  we  write 

|  u(x,y)q(x,y)dxdy  -  [u,Bq] 
wing 

The  right  side  is  considered  as  a  scalar  product.  The  matrix  B 
arises  in  the  integration  process  (loosely  speaking,  it  is  the 
substitute  for  the  area  element). 


With  this  notation  we  rewrite  Eq .  (59)  cast  in  the  form 
(60).  The  potential  41  is  expressed  in  terms  of  h  by  Eq . 

/s 

(57)  and  h  is  expressed  by  Eq .  (56);  h'(0),  which  arises  when  one 
forms  h'(t)  from  Eq.  (56)  by  differentiating  with  respect  to  the 
upper  limit  of  the  integral  is  found  in  Eq .  (58b).  Then 

j  t  A.. 

||  Pa  (t  ,x  ,y  )qn(x  ,y  )dxdy  -  +2tt{M  (c  l  ( t )  [qn  ,  Bw^  ]  +  |  c^  (  0  )  [qn  ,  Bh  £  ( t-o  )  ]d  0  ) 
wing  o 

*  /  "trailing’*,  (xtrai  llng^  >  ^ 
trailing  edge 

t  +  ~ 

-  j  °jl(o)[qn  xBhil(t-a)do)  ]do} 
o 

A  A 

The  expression  [qn,Bh^*]  and  can  be  evaluated  in 

conjunction  with  the  evaluation  of  h^ .  If  one  does  not  need  to 

examine  the  pressure  distribution  by  itself,  this  will  reduce  the 
storage  requirements  considerably. 


m 


SECTION  XII 

TRANSITION  TO  A  SMALLER  VECTOR 

Because  of  the  size  of  the  matrix  A  (in  Eq .  39),  the 
integration  of  the  integro-dif f erential  equation  for  h  is 
probably  best  carried  out  by  a  predictor-corrector  method,  rather 
than  by  some  implicit  procedure.  (The  vector  h  may  have  several 
hundred  components.)  An  explicit  method,  like  a  predictor 
corrector  method,  becomes  unstable  if  the  time  step  is  too  large. 
The  limit  of  the  time  step  which  arises  in  this  manner  is 
dependent  upon  the  discretization  in  space.  Incidentally, 
according  to  the  experiences  gained  for  partial  differential 
equations  without  memory  (as  contrasted  to  the  present  problem 
where  memory  terms  are  introduced  by  the  retardation) ,  the 
approximation  to  the  original  problem  (that  is  the  problem  before 
the  discretization)  is  poorer  if  one  uses  a  time  step  which  is 
smaller  than  the  maximum  time  step  compatible  with  stability.  In 
the  first  stages  of  the  response  to  a  function  w  whose  time 
dependence  is  given  by  a  step  function,  this  restriction  to  a 
small  time  step  is  no  disadvantage  because  a  large  time  step 
which  is  possible  in  implict  methods  will  average  out  the  rapid 
changes  of  the  potential  characteristic  of  these  initial  stages. 

After  this  initial  phase,  the  function  h(t)  changes  more 
slowly  and  finally  approaches  the  values  pertaining  to  the  steady 
state  generated  by  the  function  w^(x,y)  (here  w^) .  It  is  then 
possible  to  proceed  in  larger  time  steps  provided  that  h  is 
redefined . 

Beside  the  required  time  resolution,  the  grid  size  is 
determined  by  the  accuracy  desired  for  the  integrations  over  5 
and  n.  The  original  size  of  the  grid  will  probably  be  determined 
by  the  accuracy  desired  in  the  initial  time  response;  the  grid  so 
obtained  is  likely  to  be  sufficient  for  the  integration  with 
respect  to  E,  and  n.  The  results  obtained  for  this  fine  grid  are 
the  basis  for  the  transition  to  a  representation  of  h  by  a 
reduced  number  '  parameter^  . 

1  00 
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In  this  process,  one  can  distinguish  between  three  steps. 

In  the  first  step  one  represents  h  by  a  reduced  number  of 

parameters.  This  leaves  the  number  of  points  at  which  the  upwash 

is  evaluated  unchanged,  and  with  it  the  number  of  conditions  that 

■> 

can  be  imposed  to  this  form  of  h.  In  the  second  step,  one 
reduces  the  number  of  conditions  derived  from  the  integral 
equation,  so  that  it  matches  the  number  of  parameters  used  to 
characterize  h.  In  the  third  step,  one  enlarges  the  time 
interval. 


We  assume  immediately  that  the  time  dependence  of  w  is  given 
by  a  step  function.  Then  the  function  h  approaches  a  steady 
state  and  changes  more  and  more  slowly  after  an  initial  phase. 

We  consider  the  equation  obtained  after  a  discretization  in  space 
but  without  discretization  in  time 


■|^  +  |  A(  t  )h  ( t-T  )dT  =  w 
o 


where  w  is  constant  for  t  positive.  Assume  that  the  initial 
phase  of  the  time  evolution  of  h  during  which  this  equation  is 
used  in  its  original  form  terminates  at  time  t^.  Then  w r  write 


t-t. 


1  t 

^  +  j  A(t ) h ( t- 1  )dT  +  j  A( t )h ( t-t )di  =  w 

o  t-t1 


In  the  second  integral,  the  argument  of  h  varies  from  t^  at  the 
lower  limit  to  0  at  the  upper  limit.  For  these  values  of  t,  h  is 
already  known.  Presumably,  it  changes  fairly  rapidly  in  space  as 
well  as  in  time.  Therefore,  one  cannot  economize  in  the 
evaluation  of  this  integral.  One  will  remember,  however,  that 
the  elements  of  the  matrix  A  are  zero  if  the  retardation  t 
exceeds  a  certain  value  (which  depends  upon  the  dimensions  of  the 
wing).  Exceptions  are  matrix  elements  which  give  the  influence 
of  the  trailing  edge  on  the  upwash  at  other  parts  of  the  wing, 


because  the  trailing  edge  data  incorporates  information  about  the 
wake,  which  theoretically,  at  least,  extends  to  infinity. 

The  transition  to  fewer  parameters  in  the  representation  of 
h  can,  however,  be  made  in  the  first  integral  of  the  last 
equation.  With  suitably  chosen  functions  f£^(x,y),  one  sets 

*1 

Ml 

h ( t ,  x , y )  =  I  a. ( t  )f  J  (x,y) 

2,  =  1  1 

By  making  the  transition  from  h(t,x,y)  to  h(t),  one  obtains 

h(t)  =  M^a(t) 


where  the  matrix  M  ^  has  as  column  vectors  the  functions  f^1^ 
evaluated  at  the  control  points 


,(D 


f^Cx-.y.) 


fi1)(xi*yi) 


and 


(61  ) 


a(t )  = 


a1  (t) 
a'(t) 


Here  i  is  the  subscript  for  a  control  point  in  a  one-dimensional 
numbering  system  (and  x.,y.  are  its  coordinates).  One  thus 
arrives  at  the  equation 


t-t. 


M 


( 1  )  d  a 
dt 


[  A ( -t  )M  ^ 1  ^  ]a  ( t- t  ) dt  ) 


+  /  A( t )h ( t-t  )di 
t-t. 


=  w 


(62) 


With  this  equation,  one  can  evaluate  the  upwash  at  all  control 
points.  The  naj  dependent  "ariable  is  the  vector  a.  The  vector 


a  is  considerably  smaller  than  the  vector  h.  Of  course,  this 
equation  cannot  be  satisfied  exactly  at  all  control  points. 

In  a  second  step,  one  premultiplies  the  last  equation  by  a 
( 2 ) 

matrix  M  with  JL ^  rows  and  a  number  of  columns  which  matches 
the  dimension  of  h.  One  then  obtains 


t-t. 


M 


(3)§f  .  | 


)an ( t-x  )dx 


(63) 


t 

.  1 
t-t. 


M(5)  h  (t-T  )dx  =  w(1) 


where 


m(3)  = 


■(<*) 


'(x) 


M(2) A(x )M( 1 } 


M 


(5) 


w 


(x  ) 
(1 ) 


M(2)A(t  ) 


M(2)w 


,(2) 


The  matrix  M  may  single  out  a  suitable  number  of  control 
points,  or  its  rows  may  be  given  by  the  values  of  weight 
functions  at  the  control  points.  In  the  last  form,  the  equation 
holds  for  any  value  of  t  >  t^  In  the  description  of  the  program 
we  assumed,  however,  that  one  deals  with  equidistant  values  of  t 
and  that  the  interpolations  needed  are  incorporated  into  the 


matrix  C 


( i j ) , (mn) ,k  ’ 


There  the  numbering  of  the  control  points 


by  pairs  of  subscripts  is  used;  for  the  rows  (ij),  for  the 
columns  (m,n).  Introducing  an  analogous  notation  in  the  matrices 
one  has  for  the  matrix  elements 


M( 1 )  and  M(2) 


,o ) 


(mn )  , 


.11)/ 


vxmn’ymn) 
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a.*.  f -1-*.  . 


l,  (i  j) 


mu  u  m  ipui 


nwtiwir  www-jvw wm  w  ^  gr  in  iTV«y«v  «  v  i*r?77777r*ri 


M 


(2) 


f  ^  ( x  y  ) 

*£  iXij,yij; 


Then,  one  has 
(3) 


M 


M 


l 

(4) 


■  A  '•• 

>  •! 


••  V 


V  w  (  2 )  M(l) 

=  I  M .  ,  .  .  x  M  /  . 


"i,j  *"1  (  i  J  )  (  i  J  )  *>2 


y  y(2)  w  (  2  )  „  M(1) 

2- 1  2-  2 ;  k  i,j  num  *i  (ij)L(ij)(mn);k)M  (mn),l2 


(4) 

The  Matrix  M  which  previously  depended  upon  x,  now  depends 

(  4 ) 

upon  a  parameter  k.  This  is  expressed  by  the  notation  M  (k) 
(5) 

Similarly  for  M  .  There  elements  are  given  by 

M(5)  _  v  „(2)  P 

^.(mn^k  .  ^  M«,1  (ij)L(i  j)(mn)  ,k) 

(1)  -  M  ( 2 ) 

1  =i!j  h»(iJ}  (ij) 


Eq.  (44)  now  appears  in  the  form 


+  k  k1  ~  1  ^  (  4 )  -►  k 

(3)(da)  ♦  I  mJ  ;a  ♦  E 

n  ^arr;k  «. -o  %  ^  ?,=k-k 


M(1  }h 

hk-  £ 


1 


w(1);  k  >  k. 


d  3 

The  derivative  (-tt-V  is  obtained  by  premultiplication  with 
M ( 3 )  - 1  .  dt  k 


t  da, 

Wk  ♦ 


k'Er,)M(6)a  *  E  MmS 

£=0  1  k ~l  £=k-k1  1  k_i 


w(2) 


wi  th 


m£6)  - 

.  u(3)-l„<5) 


.t -*l.  ..  .1 


.  aL  iw  ii.ii.  iw  il.  il.il.il.  iU  al.il.  it.  iU  il.  it.  tWiLiLiU  iLil.  iL  iLiU  ik  fL'Ak.il.'tLil.'iLi 


i5(2) 


m(3)-^(D 


In  addition,  one  needs  initial  conditions  for  the  vector  a  at  k  = 
k,  .  Known  from  the  preceding  computation  is  h  for  k  =  k...  One 
must  approximate  h(k^)  by  M  a(k^) 

This  is  the  formulation  for  the  original  time  step.  A 
reduction  in  the  number  of  parameters  describing  h(t,x,y)  is 
likely  to  tolerate  a  larger  time  step.  Such  a  transition  can  be 
made  if  the  new  time  step  is  a  multiple  of  the  original  one.  One 
computes  the  vector  a  only  for  values  belonging  to  the  larger 
time  step  and  expresses  intermediate  values  of  the  vector  a  by 
linear  interpolation.  We  show  the  procedure  for  tripling  the 
time  step.  (It  might  be  more  realistic  only  to  double  the  time 
step.  Tripling  is  chosen  to  indicate  how  one  proceeds  in 
general.)  Then  only  values  for  k  -  k1  +  3s ; s  =0,  1...  occur. 

Let 


k1  +3s 


In  the  second  sum  of  the  last  equation,  the  summation  extends 


for  s  »  0,  (k«k.|),  from  l  •  0  to  l  ■  k1 


for  s-1,  (k-k.j+3),  from  2,  ■=  3  to  l  «  kj+3 


for  s  =  2,  (k  =  ^+6)  from  l  *  6  to  l  >  k.j+6 

In  all  cases  the  subscript  of  k-t  of  h  varies  from  0  to  k^.  The 

number  of  terms  in  the  sum  is  always  the  same  and  does  not  depend 

upon  the  size  of  the  new  integral.  However,  for  s  sufficiently 

( 5 ) 

large,  all  elements  of  the  matrix  M  vanish  except  for  those 
pertaining  to  values  of  h  at  the  trailing  edge. 

In  the  first  sum,  the  interpolation  for  the  vectors  ak_ ^  is 
introduced.  One  ha l 


A  A  A  »\ 


V.3S.1  '  <2/3>as  *  l1/3)as.1 


V3..2  •  (,/3)as  *  <2/3)as.1 


The  summand  of  the  first  sum  are  now  ordered  by  the  subscript  of 
-> 

/V 

a.  For  (k  =  ^+3),  one  obtains  three  terms  (fc-O,  1-1,  l«2) 
^Q6)a.,  +  M1(6)((2/3)a1  ♦  (1/3)aQ)  ♦  M^6)((1/3)a1  +  (2/3)aQ) 

-  (M<6)-  (2/3)M1(6)  +  (1/3)M^6)  )a1  ♦  ((1/3)M1(6)>-(2/3)M26)  )aQ 
For  s  -  2  (k  -  k^+  6),  one  obtains 
mJ6)+  ( 2/3)m| 6) +  (1/3)M^6))a2 

♦  ((1/3)M1(6)+(2/3)M^6)+  M^6)  +  (2/3)mJ6)+  (1/3)M^6))a1 

♦  ((1/3)mJ6)  ♦  (2/3)M^6))aQ 

For  general  s  (k  »  k..  ♦  3„),  one  obtains 

I  O 


^6>a3  *  «i6)aV,  * 


i‘6)  a  ...Ht6**>;n 
n  s-n  0 


with 


mq6)  =  +  (2/3)M1(6)+  (1/3  M<6) 


l^6)  -  (1/3)M^)0  *(2/3)m15)-  ♦  M^)  +  (2/3)M^?1  ♦  (1/3)M~ 


These  quantities  are  independent  of  s  in  contrast  to 


M(6,s)  -  (1/3)M^*2  +(2/3)M^l1 


The  first  sum  then  becomes 


#i’  n  *  »i6,3>^ 

n-0  n  s-n  3  0 


and  one  obtains 


.  V  «<6>a,  „  .  A«^>L 

dt  n_Q  n  s-n  0 


k0+3s_1  (7) 

+  j, » 3s  M  1  V^s-t  "  W ( 2 ) 


(65) 


The  matrices  M  and  the  vectors  and  w^2^  can  be  evaluated  in 

advance.  For  s  =  0,  the  first  sum  is  zero. 
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SECTION  XIII 

EXAMPLES  FOR  THE  TRANSITION  TO  A  SMALLER  VECTOR 


The  procedure  is  illustrated  by  two  examples  for  the  two- 
dimensional  case.  With  some  modifications  the  ideas  can  be 
carried  over  to  the  three-dimensional  case. 

We  assume  that  control  points  have  been  chosen  along  the 
chord  of  the  profile  at  10  equal  intervals.  (One  may  find  it 
desirable  to  use  more  control  points,  but  10  is  sufficient  to 
illustrate  the  procedure.)  In  the  first  example,  the  size  of  the 
interval  is  doubled,  the  amount  of  labor  to  proceed  over  an  equal 
time  span  is  then  reduced  by  a  factor  of  1/4. 

The  function  h  is  now  expressed  by  its  value  of  only  five 
control  points.  For  the  value  of  h  at  the  control  points  that 
have  been  left  out  are  determined  by  linear  interpolation  between 
the  neighbors.  At  the  leading  edge,  one  has  a  square  roots 
singularity.  This  is  taken  into  account  when  h^h  at  control 
point  1)  is  expressed  by  the  value  of  h2-  One  then  has 


1/2  1/2 


1/2  1/2 


The  condition  for  the  vector  a  are  imposed  at  the  new  control 
points.  Accordingly 


0  0 


0  0 


0  0 


0  0 


Then 


)  «  M(3) 


where  1^  is  the  5  by  5  identify  matrix 

^  4-  U  /.a  _  U  (  ^  ) 


One  then  evaluates  M 


M(6)  and  M(5)-  M(7) 


For  the 


two-dimensional  case,  the  rows  and  columns  of  the  matrices  C  are 
numbered  by  single  subscripts. 


The  integration  with  respect  to  t  can  now  be  carried  out 
with  twice  the  original  time  interval 


k  -  k1  +  2s 


Then,  one  has  r.  following  formulae 


F\«w»wnj»uwwwwww  wiwt  imnnr*  p»w*twiruinw^y  wn^jr*ywjcmjrmjcmjrrjrw7fj,j^frfjv£fu A.  wv ^ ~-r j 


m!6)  -  Mq6)  +  (1/2)M® 


M 


0 

(6) 


<1'2>M 2n-1  *  M2n’  *  (,/2)H2n!l 


M(6,S)  -  <1/2)M2s-1 


The  initial  conditions  for  a  are  the  value  of  h(k.j)  at  the  new 
control  points 


^ an ^o  “  (h2n)k. 


The  problem  then  appears  in  the  form  as  shown  in  Eq .  (65) 


J 


•1 


S-  1 

S2s  .  (  m' 

dt  J  n 

n  =  0 


<6>;  .  .  *°T~'  h”>  h„  ,  , 

3-n  o  t.2s  1  k0.2a-l 


w 


(2) 


Further  details  are  left  the  the  reader. 

This  procedure  i3  somewhat  crude  insofar  as  it  expresses  h 
at  the  intermediate  control  point  by  linear  interpolation  between 
the  control  points  except  for  the  interval  next  to  the  leading 
edge.  While  the  approximation  for  h  by  a  piecewise  linear 
function  in  10  intervals  may  be  sufficient,  it  may  not  be 
sufficiently  accurate  if  one  has  only  5  intervals.  The  approach 
is  sensible  if  the  original  interval  had  been  chosen  because  of 
the  desire  to  obtain  a  close  time  resolution  while  it  is 
unnecessarily  fine  for  the  representation  of  the  final  steady 
state.  In  the  two-dimensional  case,  the  labor  required  to 
proceed  by  the  same  time  is  reduced  by  a  factor  of  1/4;  in  the 
three-dimensional  case  by  a  factor  of  1/8.  One  can,  of  course, 
introduce  a  matrix  in  which  one  uses  higher  order 

interpolation  formulae  for  h.  The  next  procedure  is  of  this 
character  insofar  as  uses  for  the  f^1^  analytical  expressions 
applied  over  the  whole  profile. 

For  this  purpose,  one  can  use  functions  derived  from 
linearized  airfoil  theory.  In  the  three-dimensional  case,  one 
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may  use  such  expression  for  the  chordwise  dependence  of  h,  their 
coefficients  are  then  functions  of  a  spanwise  coordinate. 

Classical  two-dimensional  airfoil  theory  is  partially 
summarized  in  Appendix  E  of  Reference  1.  The  profile  is 
represented  by  slit  in  the  x,y  plane  extending  along  the  x  axis 
from  -1  to  +1.  The  x,y  plane  with  this  slit  is  mapped 
conformally  into  the  exterior  of  a  circle  with  radius  1.  At  this 
circle,  the  potential  (in  this  case  the  function  h), 
antisymmetric  with  respect  to  the  x  axis,  is  expressed  by  a 
complete  system  of  functions 


v 

•j 

iA 


hi  =  -sin  6 
n 


.-I 


where  9  is  the  polar  angle.  For  the  mapping  from  the  circle  to 
the  slit  in  the  x,y  plane  one  has 


x  =  cos  9  ;  (y=0) 

Thus,  a  system  of  functions  suitable  to  represent  h  (without 
allowing  for  circulation)  is  given  by 

hn  =  -sin  (n8) 


r 


with 

sin  9  -  ( 1  -x2  )1  /2 

To  these  expressions  of  h,  one  adds  a  circulation  term  so  that 
for  0=0,  which  is  the  map  of  the  trailing  edge  dh/d0  =  0.  (One 
remembers  that  dh/de  in  the  circle  plane  is  the  tangential 
derivative  of  h  along  the  circle.)  This  condition  ensures  that 
in  the  x,y  plane  the  velocity  at  the  trailing  edge  is  finite. 

Thus 


h  =  -  sin(ne)  +  n(  0-tt  ) 
n 

1  1  1 
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The  upwash  pertaining  to  such  a  function  hn  at  the  slit  in  the 
x,y  plane  is  found  according  to  Appendix  E  of  Reference  1  as 


wn  =  n  sin(n8)/sin  0 


The  upwash  generated  by  the  circulation  term  is  zero.  The 

expression  wn  is  a  polynomial  in  x.  If  the  function  w  prescribed 

in  the  program  for  the  unsteady  flow  is  such  a  function  w  and 

n 

its  time  dependence  is  given  by  a  step  function,  then  the 
asymptotic  response  wi! 
orthogonality  relation 


asymptotic  response  will  be  the  function  h  .  One  has  the 


M 

/  sin(ne)  sin(n0)d6  =0;  m  ^  n 


hence 


But 


M 

I  w„  w 


n  m 


sin  edo  =0,  m  4  n 


a 

A 

A 

>, 

'A 

A 


Hence 


+  1 


cos  6  =  x 


-  sin  edo  -  dx 


\  wm  (1-x2)1 /2dx  -  0,  m  4  n 


-1 


n  m 


By  this  formula  one  can  decompose  an  upwash  given  by 


n 

w  .  I  a  w„ 
^  n  n 
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into  the  individual  components,  which  are  generated  by  the 

functions  h  . 

n 

During  the  transition  to  the  asymptotic  value,  the  wake  is 
not  fully  developed.  Therefore,  the  functions  hn  by  themselves 
are  not  quite  sufficient  to  approximate  the  potential.  The 
individual  expressions  give  at  the  trailing  edge 

hn(x)  -  const  +  (x-xt^il 

therefore,  h  =0.  But  at  the  trailing  edge,  one  must 
n  |  x 

prescribe  the  conditions  that  no  pressure  jump  occurs  between  the 
upper  and  lower  side  of  the  profile.  Therefore, 

*x  +  M  1 4>t  =  0 

for  the  trailing  edge.  To  allow  for  this  expression,  we  include 
in  the  list  of  functions  used  for  the  representation  of  h,  a 
f unot i on 

WN+1  *  (x+1 ) 1 72 ( x~ 1 ) 

At  the  trailing  edge  it  gives  a  zero  contribution  to  h  and  a 
1/2 

value  of  2  for  the  value  of  hx.  At  the  leading  edge,  it  has 
as  all  other  functions  a  square  root  singularity.  This  is  then 
the  set  of  functions  f^  used  to  express  h  in  terms  of  the 
parameters  an. 

( 2 ) 

A  matrix  M  chosen  in  such  a  manner  that  in  the  steady 
case  it  gives  the  above  decomposition  of  the  downwash,  is  likely 
to  be  pratical.  Therefore 

f^2^  -  wn(1-x2)1 72  -  sin(ne);  n*1 . N 

( 2 ) 

This  defines  N  of  the  functions  f  which  appear  in  the  row 

(  2 ) 

vectors  of  the  matrix  Mv  .  One  needs  an  additional  condition, 
for  this  we  choose  the  requirement  that  the  original  discretized 


expression  for  the  upwash  is  satisfied  at  the  trailing  edge.  The 
trailing  edge  is  singled  out,  because  of  its  unique  role  in 
determining  the  vortex  shedding  in  the  wake. 

To  carry  out  the  procedure  for  the  new  representation  for 
h,  one  must  have  the  initial  conditions  for  the  vector  a^.  It 
is  assumed  that  the  procedure  in  its  original  form  has  been 
carried  out  up  to  a  station  k1  in  time  so  that  the  vector  hk1  is 
known.  Then  one  must  find  the  values  of  an  which  give  an 
acceptable  approximation  for  h^ 

h.  -  £  a  h  ♦  aM  ( x 1  )  '^(x-1  ) 
k .  ^  n  ,  o  m  N  + 1 

From  the  known  values  of  hk1 ,  one  determines  numerically  hx 
at  the  trailing  edge.  Then,  because  hm  x  =  0 

J/2 

aN+1,0  c  nx  trailing  edge 
Now  we  are  left  with 


N 

£ 

1 


where  one  writes  the  last  expression  as  a  vector  with  components 
obtained  by  inserting  the  values  of  x  for  the  original  points. 


There  are  orthogonality  relations  for  the  functions  h  if 

w 

one  leaves  out  the  contribution  of  the  circulation.  The  combined 
contribution  of  the  circulation  at  the  trailing  edge  gives  hk1  at 
the  trailing  edge.  Therefore 


N 

£ 

1 


a  h 


no  n 


WxM) 


(x-1  ) 


’kl  ,  trail  ing  edge 


where  the  values  of  8  are  re-evaluated  at  the  control  points. 


The  orthogonality  conditions  for  the  hn  are 


/  sin  6  sin  (ra  0)d6  =  0 


2  1/? 

/  sin  9  sin  (m  6)  (1-x  )  dx 


sin  (m&)/sin  0  is  a  polynomial  in  x.  Since  hk1  is  given  only  at 

the  control  points,  one  must  use  some  simple  approximate 

integration  formula  if  one  applies  the  orthogonality  condition. 

The  matrix  then  obtained  for  computing  the  coefficients  anQ  is 

nearly  diagonal.  Therefore,  it  is  rather  simple  to  obtain  the 

values  of  a  .  If  one  uses  these  values  of  a  and  evaluates  the 
no  no 

potential  at  the  trailing  edge  given  by 


E  a„  nn 
1  no 


one  will,  in  general,  not  obtain  the  original  value  of  hk  ^  at  the 

trailing  edge.  It  may  be  desirable  to  have  continuity  for  the 

potential  at  the  trailing  edge  between  the  values  computed  by  the 
original  procedure  and  that  computed  with  the  vectors  a.  Then 
one  must  take  some  correction  to  the  values  of  a 

no 

With  this  procedure,  it  may  be  possible  to  carry  out  a  more 

radical  reduction  in  the  size  of  the  problem.  And,  it  may  also 

be  possible  to  apply  a  much  larger  time  step. 

Other  points  of  view  than  those  anticipated  here,  may,  of 
course,  emerge  in  the  course  of  numerical  experimentation. 
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SECTION  XIV 

THE  METHOD  IN  THE  FRAMEWORK  OF  AN  AEROELASTIC  PROBLEM 


The  motivation  for  the  development  of  the  method  described 
here  has  been  the  treatment  of  unsteady  aeroelastic  problems  in 
the  physical  space  rather  than  in  the  frequency  space.  This 
determines  the  manner  in  which  the  results  will  be  used  and  the 
required  accuracy.  We  discuss  this  in  some  detail. 


The  basic  equations  of  aeroel as ti ci ty  express  the 
equilibrium  between  the  elastic  forces,  the  damping  forces  (both 
within  the  structure),  the  acceleration  forces  (d'Alembert 
forces),  and  the  aeroelastic  forces.  All  these  forces  are 
expressed  in  terms  of  the  displacements  of  the  wing,  in  this 
report  denoted  by  g(t,x,y).  The  elastic  forces  are  expressed  in 
terms  of  the  g  itself,  the  damping  forces  by  3g/3t,  the 


2  2 

acceleration  forces  by  3  g/3t  ,  all  expressions  taken  at  the 


current  time,  the  aerodynamic  forces  are  expressed  by  3g/3t  and  g 
at  the  current  and  at  retarded  times.  One  can  consider  g(t,x,y) 
at  fixed  time  t  as  the  element  of  a  function  space,  called 


2  2 

displacement  space;  3g/3t  and  3  g/3t  are  then  other  elements  of 


this  space.  If  t  varies  these  elements  will,  of  course,  move 
along  curves  in  the  displacement  space.  The  elastic,  damping, 
mass,  and  aerodynamic  forces  are  obtained  from  these  elements  by 
operators,  which  map  these  elements  from  the  displacement  space 
into  a  "force"  space. 


The  equations  of  motions  express  the  requirement  that  the 
sum  of  the  forces  generated  by  these  mappings  give  at  all  times 
the  zero  element  of  the  force  space. 


If 
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restricts  in  the  process 
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Accordingly  the  sum  of  these  forces  will  not  vanish  perfectly. 

An  exception  occurs,  if  one  restricts  oneself  to  the  mass  and 
stiffness  operators.  The  mass  and  stiffness  operator  operat''ng 
on  functions  of  x  and  y  define  an  eigenvalue  problem.  For  the 
eigen solutions  the  mapping  due  to  the  mass  and  the  stiffness 
operator  give,  except  for  a  factor  of  proportionality  (the 
eigenvalue),  the  same  vector  in  force  space.  If  the  subspace 
used  for  the  approximation  of  the  deformations  is  spanned  by  a 
finite  number  of  eigenvectors  then  the  subspaces  obtained  by 
applying  the  mass  and  stiffness  operator  are  identical.  By 
defining  the  spaces  in  a  somewhat  different  manner  and  admitting 
complex  elements  one  can  also  include  the  damping  operator.  The 
aerodynamic  operator,  however,  is  different  for  each  value  of  the 
retardation  and  cannot  be  included  in  such  an  approach.  If  one 
had  only  a  stiffness  and  a  mass  operator  and  used  a  subspace 
spanned  by  eigenfunctions,  then  one  would  satisfy  the 
differential  equations  of  the  homogeneous  problem  perfectly  (an 
approximation  arises  only  in  the  manner  in  which  one  satisfies 
the  initial  conditions). 

As  mentioned  before,  the  aerodynamic  forces  are  c.ctermined 
by  operators  which  are  different  for  each  value  of  the 
retardation,  and  therefore  do  not  lend  themselves  to  such  an 
approach.  If  one  applies  a  discretization  in  displacement  space, 
then  one  cannot  satisfy  the  differential  equations  of  motions 
(which  are  relations  in  force  space)  perfectly. 

The  problem  is  ultimately  reduced  to  a  system  of  ordinary 
differential  equations  for  the  coefficients  of  the  spanning 
vectors  of  v.he  deformation  space.  To  obtain  such  a  system,  one 
first  expresses  the  displacement  vector  (the  elements  of  the 
displacement  subspace)  by  the  sum  of  spanning  vectors  of  the 
displacement  space  multiplied  by  time  depend  coefficient  (so  far 
unknown)  and  substitutes  them  into  the  equations  of  motion  (in 
other  words  one  applies  to  them  the  stiffness,  damping,  mass,  and 
aerodynamic  om  ’atoru ) .  This  gives  a  vector  in  the  force  space. 
Subsequently,  one  muitipliis  the  equations  of  motion  by  weight 
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functions  and  then  integrates  over  x  and  y.  This  generates  a 
system  ordinary  differential  equations  for  the  coefficients  of 
the  spanning  vectors.  In  problems  where  the  concept  of 
eigenvectors  is  applicable  there  is  no  question  about  the  choice 
of  the  weight  functions,  they  are  the  spanning  vectors  of  the 
displacement  space.  Because  of  the  orthogonally  relations  which 
exist  in  this  case,  one  even  obtains  single  equations  for  the 
individual  coefficients.  The  simple  example  of  an  oscillating 
string  with  non-uniform  mass  distribution  can  serve  as  an 
illustration  for  this  situation. 

Such  thorough  simplifications  are  not  feasible  in  the 
presence  of  aerodynamic  forces.  The  eigenfunctions  of  the 
problem  without  aerodynamic  forces  lose  their  physical  and 
mathematical  significance,  although  one  can  still  expect  that  the 
eigenfunctions  pertaining  to  the  lowest  frequencies  are  a 
"reasonable"  spanning  of  the  displacement  space.  Instead  of 
determining  these  eigenfunctions,  one  describes  the  displacements 
by  a  finite  number  of  linearly  independent  reasonably  smooth 
functions.  If  one  chooses,  for  instance  10  su.  h  functions,  then 
it  is  quite  likely  that  a  linear  combination  of  them  will  give 
good  approximations  to  the  first  few  eigenfunctions,  and  this  is 
sufficient.  If,  in  a  specific  situation,  the  coefficients 
ultimately  obtained  for  those  functions  which  have  the  greatest 
waviness  are  small,  then  the  choice  of  the  system  of  spanning 
functions  is  "reasonable." 

The  weight  functions  applied  to  the  force  equations  in 
order  to  arrive  at  a  system  of  ordinary  differential  equations 
for  the  coefficients  are  chosen  in  analogy  to  the  treatment  in 
terms  of  eigenfunctions,  i.e.,  one  chooses  the  functions  which 
span  the  displacement  subspace.  By  the  subsequent  integration 
over  x  and  y  one  obtains  the  stiffness,  damping,  mass  and  the 
aerodynamic  matrices.  The  latter  depend  upon  the  retardations. 
The  others  are  time  independent.  The  weight  functions  applied 
here  are  the  same  as  the  functions  q  mentioned  in  Section  XI.  If 
the  stiffness  damping  and  mass  matrix  are  determined 
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experimentally,  then  one  must,  of  course,  use  in  the  processing 
of  the  aerodynamic  forces  the  weight  functions  of  the 
experimental  arrangement.  There  is  nothing  new  in  this 
description,  every  flutter  computation  is  based  on  this  idea. 

The  model  obtained  after  one  has  introduced  finite 
stiffness  damping  and  mass  matrices  defines  certain 
ei genf requenci es  and  by  them  a  measure  for  the  time  resolution 
which  is  of  technical  interest.  This,  in  turn,  defines  the 
maximum  grid  size  by  which  one  achieves  the  desired  time 
resolution.  (Another  criterion  for  the  upper  bound  of  the  grid 
size  is  the  accuracy  of  the  integrations  which  appear  in  the 
determination  of  the  aerodynamic  force.) 

To  illustrate  the  relation  between  time  interval  and 
frequency  we  consider 


a”  +  v2a  =  f(t) 

This  equation  can  be  interpreted  as  the  equation  for  the 
coefficient  of  one  eigenfunction  under  the  influence  or  a  forcing 
function  f(t);  which  takes  the  place  of  the  aerodynamic  forces. 

The  solution  (with  initial  conditions  a  =  0,  a’  0)  is 
given  by 

_1  t  t 

a(t)  =  v  {sin  ( vt)/f (t )cos( vt )dt  -  cos (vt ) / f (t ) si n ( vt )dt } 

o  o 

If  f(t)  is  of  significant  size  only  in  the  internal  0  <  t  <  t, 

- 1  1 
t1 <<v  ,  then  one  obtains  as  approximation  for  t  >  t1 

- 1  ^ 

a(t)~v  sin(vt)  /f(t)dt 

o 

In  other  words,  it  does  not  matter,  if  in  an  interval  At<<v  1  the 
function  f(t)  is  known  only  in  the  average.  One  might  consider 
At  =  v"1  10-1  as  sufficient. 
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The  aerodynamic  forces  are  governed  by  two  phenomena  which 
may  have  different  time  scales.  One  is  the  propagations  of  waves 
over  the  wing.  We  have  considered  a  function  w(t,x,y)  whose  time 
dependence  is  given  by  a  step  function.  At  the  time  immediately 
after  the  step,  the  pressure  at  the  wing  is  solely  determined  by 
the  local  value  of  w(x,y).  But  this  pressure  distribution  is  not 
compatible  with  the  pressures  outside  of  the  wing  (which  are 
zero).  The  pressure  adaptation  therefore  needed,  occurs  by  waves 
which  travel  over  the  wing.  In  the  two  dimensional  case  one  has 
only  waves  traveling  in  upstream  and  downstream  directions.  The 
downstream  travel  occurs  at  higher  velocities  than  the  upstream 
travel,  especially  at  high  subsonic  Mach  numbers.  In  the  three- 
dimensional  case,  one  has  waves  in  all  directions.  These  waves 
are  reflected  at  the  wing  contour,  but  they  attenuate  with 
distance  because  the  underlying  phenomena  space  are  three 
dimensi onal . 


A  second  time  scale  is  given  by  the  development  of  the 
wake.  While  the  adaptation  to  the  outside  pressures  is 
practically  over  after  a  certain  (possibly  rather  short)  time, 
the  wake  keeps  extending  and  changing.  Especially  at  low  Mach 
numbers  the  development  of  a  steady  wake  will  take  a  considerable 
time.  The  wake  vorticity  affects  the  upwash  on  the  wing  and  with 
it  the  pressure  distribution.  Some  time  after  the  step  in  the 
upwash,  the  upwash  at  the  wing  generated  by  the  wake  will  change 
only  slowly  and  the  adaptation  to  this  upwash  will  probably  occur 
in  a  quasisteady  manner.  In  testing  the  method  one  should  be 
aware  of  this  distinction  and  try  to  recognize  the  relative 
importance  of  these  phenomena.  The  characteristic  times  for 
these  phenomena  should  be  considered  in  relation  to  the  time 
resolution  imposed  by  mechanical  properties  of  the  system.  The 
transition  from  a  finer  to  a  coarser  grid  described  in  the 
preceding  section  is  the  practical  consequence  of  the  difference 
between  the  initial  phase  where  quickly  moving  waves  determine 
the  pressure  distribution  and  a  later  phase  where  such  waves  play 
a  minor  role  and  the  phenomena  are  quasisteady. 
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SECTION  XV 


EFFECT  OF  A  TRUNCATION  OF  THE  WAKE 


The  wake  generated  by  an  upwash  w(t,x,y)  whose  time 
dependence  is  given  by  a  step  function,  starts  with  length  zero. 
The  end  of  the  wake  moves  downstream  with  the  velocity  of  the 
particles,  that  is  with  the  free  stream  velocity.  If  one  follows 
the  evolution  of  the  flow  field  for  a  long  time  then  the  wake 
will  become  very  long  and  the  labor  to  take  the  distant  parts 
into  account  will  be  excessive,  while  their  influence  on  the 
upwash  at  the  wing  will  become  smaller  and  smaller.  In  practice 
the  wake  will  be  terminated  after  a  finite  distance. 

We  try  here  to  describe  the  effect  of  such  a  truncation  in 
general  terms.  Actually  this  is  done  only  for  the  final  steady 
si  '  e ,  which  arises  in  response  to  a  step  function  in  time,  for 
the  interpretation  will  be  based  on  the  concept  of  vortices  in  a 
steady  flow. 

i" or  the  two-dimensional  case  h  is  zero  by  definition.  In 

y  r  2 )  ( ■)  \ 

the  functions  and  Bg  one  encounters  factors  h  and  h  (in 

a  in  )rc  familiar  notation  h  and  h  ).  Disregarding  distant  parts 

x  y 

o;  the  wake  means  that  one  sets  h  and  h  equal  to  zero.  In  a 

x  y 

wake  there  is  no  pressure  difference  between  the  upper  and  lower 
sides  in  a  steady  linearized  flow  <f>x  and  with  it  hx  are  zero.  By 
truncating  the  wake,  one  sets  hx  equal  to  zero  prematurely,  but 
in  the  final  stage  hx  is  zero  in  any  case.  Therefore,  the  final 
steady  state  for  a  two-dimensional  flow  remain  unaffected  by  the 
truncation  of  the  wake. 

In  the  three-dimensional  steady  case  one  has  vortices 
trailing  along  the  stream  lines  (here  lines  in  the  direction  to 
the  x  axis)  out  to  infinity.  Therefore,  hy  /  0  even  in  the 
steady  case.  h  is  the  local  vorticity  strength  in  the  wake,  the 

y  + 

vorticity  vector  has  the  x  direction.  But  every  vector  field  u 
satisfies  div  curl  u  =  0.  Vorticity  lines  cannot  end  somewhere 
in  the  interior  of  a  vector  field.  It  seems  as  if  this 


requirement  is  violated,  if  one  disregards  the  values  of  <J>  in 
some  parts  of  the  flow  field  further  downstream. 

The  velocity  field  due  to  a  vortex  can  be  evaluated  by  means 
of  the  Biot-Savart  law.  In  this  manner  one  can  evaluate,  for 
instance,  the  velocity  field  of  a  vortex  extending  along  the  x 
axis  from  negative  to  positive  infinity;  in  this  manner  one 
obtains  the  so-called  potential  vortex.  If  one  extends  the 
integration  from  negative  infinity  to  a  finite  value  x  =  x  ,  then 
it  seems  as  if  one  violates  the  condition  div  curl  u  =  0  at  this 
point  because  one  does  not  continue  the  vortex.  Actually  one 
creates  a  field  in  which  vortex  lines,  equally  distributed  over 
all  directions,  emerge  from  the  point  x  =  Xq  in  the  same  manner 
as  stream  lines  emerge  from  a  source.  If  one  fails  to  take  hy 
into  account  in  part  of  the  wake  for  downstream,  one  no  longer 
expresses  the  vorticity  present  in  that  part  of  the  wake  but 
introduces  vorticity  sources  at  the  cut-off  line,  which 
introduces  vorticity  in  the  entire  flow  field.  The  effect  of 
i  sc  vortices  at  a  distance  is  weak  because  the  potential 
outside  of  the  wake  is  zero.  The  integral  J  h  y d  y  evaluated  across 
the  wake  is  zero;  therefore,  one  has  for  every  vortex  source  also 
.1  vortex  sink  in  a  vicinity  smaller  than  the  vortex  width.  In 
any  case  the  final  steady  state  in  the  three-dimensional  case 
with  a  truncated  wake  is  not  the  same  as  for  a  non-truncated 
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Figure  5*  Numbering  of  the  Corners  at  and  Adjacent  to  a  Control 
Point  on  the  Axis  of  Symmetry  (j=0). 


Figure  6.  Triangle  m=0,  n,  1  (at  the  leading  edge). 
Sweep  Angle,  £  n  and  uv  Systems. 
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Figure  9.  Triangles  m=i  -1,  n,  1  and 
m=i.-l,  n,  2  Triangle  in  the 
Wake  Cn  and  uv  Systems. 
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Figure  10  Triangle  m=ifc-l,  n,  1 
uv  and  uv  Systems. 


Figure  13.  Vicinity  of  a  Control  Point  in  the 
Interior  at  the  Axis  of  Symmetry. 
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ure  1^.  Control  Point  at  the  Axis  of  Symmetry 
Next  to  the  Leading  Edge. 
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